/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Diff of /branches/diaplayground/ripley/src/Brick.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/ripley/src/Brick.cpp revision 4687 by jfenwick, Wed Feb 19 00:03:29 2014 UTC branches/diaplayground/ripley/src/Brick.cpp revision 4949 by caltinay, Mon May 19 05:54:58 2014 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
17    #include <limits>
18    
19  #include <ripley/Brick.h>  #include <ripley/Brick.h>
 #include <paso/SystemMatrix.h>  
20  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
21  #include <ripley/DefaultAssembler3D.h>  #include <ripley/DefaultAssembler3D.h>
22  #include <ripley/WaveAssembler3D.h>  #include <ripley/WaveAssembler3D.h>
23    #include <ripley/LameAssembler3D.h>
24    #include <ripley/domainhelpers.h>
25  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
26    
27  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 43  using esysUtils::FileWriter; Line 46  using esysUtils::FileWriter;
46    
47  namespace ripley {  namespace ripley {
48    
49    int indexOfMax(int a, int b, int c) {
50        if (a > b) {
51            if (c > a) {
52                return 2;
53            }
54            return 0;
55        } else if (b > c) {
56            return 1;
57        }
58        return 2;
59    }
60    
61  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,
62               double x1, double y1, double z1, int d0, int d1, int d2,               double x1, double y1, double z1, int d0, int d1, int d2,
63               const std::vector<double>& points, const std::vector<int>& tags,               const std::vector<double>& points, const std::vector<int>& tags,
64               const simap_t& tagnamestonums) :               const simap_t& tagnamestonums,
65      RipleyDomain(3)               escript::SubWorld_ptr w) :
66  {      RipleyDomain(3, w)
67    {
68        if (static_cast<long>(n0 + 1) * static_cast<long>(n1 + 1)
69                * static_cast<long>(n2 + 1) > std::numeric_limits<int>::max())
70            throw RipleyException("The number of elements has overflowed, this "
71                    "limit may be raised in future releases.");
72    
73        if (n0 <= 0 || n1 <= 0 || n2 <= 0)
74            throw RipleyException("Number of elements in each spatial dimension "
75                    "must be positive");
76    
77      // ignore subdivision parameters for serial run      // ignore subdivision parameters for serial run
78      if (m_mpiInfo->size == 1) {      if (m_mpiInfo->size == 1) {
79          d0=1;          d0=1;
# Line 56  Brick::Brick(int n0, int n1, int n2, dou Line 81  Brick::Brick(int n0, int n1, int n2, dou
81          d2=1;          d2=1;
82      }      }
83      bool warn=false;      bool warn=false;
84      // if number of subdivisions is non-positive, try to subdivide by the same  
85      // ratio as the number of elements      std::vector<int> factors;
86      if (d0<=0 && d1<=0 && d2<=0) {      int ranks = m_mpiInfo->size;
87          warn=true;      int epr[3] = {n0,n1,n2};
88          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));      int d[3] = {d0,d1,d2};
89          d0=max(1, d0);      if (d0<=0 || d1<=0 || d2<=0) {
90          d1=max(1, (int)(d0*n1/(float)n0));          for (int i = 0; i < 3; i++) {
91          d2=m_mpiInfo->size/(d0*d1);              if (d[i] < 1) {
92          if (d0*d1*d2 != m_mpiInfo->size) {                  d[i] = 1;
93              // ratios not the same so leave "smallest" side undivided and try                  continue;
94              // dividing 2 sides only              }
95              if (n0>=n1) {              epr[i] = -1; // can no longer be max
96                  if (n1>=n2) {              if (ranks % d[i] != 0) {
97                      d0=d1=0;                  throw RipleyException("Invalid number of spatial subdivisions");
98                      d2=1;              }
99                  } else {              //remove
100                      d0=d2=0;              ranks /= d[i];
101                      d1=1;          }
102                  }          factorise(factors, ranks);
103              } else {          if (factors.size() != 0) {
104                  if (n0>=n2) {              warn = true;
105                      d0=d1=0;          }
106                      d2=1;      }
107                  } else {      while (factors.size() > 0) {
108                      d0=1;          int i = indexOfMax(epr[0],epr[1],epr[2]);
109                      d1=d2=0;          int f = factors.back();
110                  }          factors.pop_back();
111              }          d[i] *= f;
112          }          epr[i] /= f;
     }  
     if (d0<=0 && d1<=0) {  
         warn=true;  
         d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));  
         d1=m_mpiInfo->size/d0;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n0>n1) {  
                 d0=0;  
                 d1=1;  
             } else {  
                 d0=1;  
                 d1=0;  
             }  
         }  
     } else if (d0<=0 && d2<=0) {  
         warn=true;  
         d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1))));  
         d2=m_mpiInfo->size/d0;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n0>n2) {  
                 d0=0;  
                 d2=1;  
             } else {  
                 d0=1;  
                 d2=0;  
             }  
         }  
     } else if (d1<=0 && d2<=0) {  
         warn=true;  
         d1=max(1, int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1))));  
         d2=m_mpiInfo->size/d1;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n1>n2) {  
                 d1=0;  
                 d2=1;  
             } else {  
                 d1=1;  
                 d2=0;  
             }  
         }  
     }  
     if (d0<=0) {  
         // d1,d2 are preset, determine d0  
         d0=m_mpiInfo->size/(d1*d2);  
     } else if (d1<=0) {  
         // d0,d2 are preset, determine d1  
         d1=m_mpiInfo->size/(d0*d2);  
     } else if (d2<=0) {  
         // d0,d1 are preset, determine d2  
         d2=m_mpiInfo->size/(d0*d1);  
113      }      }
114        d0 = d[0]; d1 = d[1]; d2 = d[2];
115    
116      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
117      // among number of ranks      // among number of ranks
118      if (d0*d1*d2 != m_mpiInfo->size)      if (d0*d1*d2 != m_mpiInfo->size){
119          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
120        }
121      if (warn) {      if (warn) {
122          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
123              << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;              << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;
# Line 230  Brick::Brick(int n0, int n1, int n2, dou Line 203  Brick::Brick(int n0, int n1, int n2, dou
203      populateSampleIds();      populateSampleIds();
204      createPattern();      createPattern();
205            
     assembler = new DefaultAssembler3D(this, m_dx, m_NX, m_NE, m_NN);  
206      for (map<string, int>::const_iterator i = tagnamestonums.begin();      for (map<string, int>::const_iterator i = tagnamestonums.begin();
207              i != tagnamestonums.end(); i++) {              i != tagnamestonums.end(); i++) {
208          setTagMap(i->first, i->second);          setTagMap(i->first, i->second);
# Line 241  Brick::Brick(int n0, int n1, int n2, dou Line 213  Brick::Brick(int n0, int n1, int n2, dou
213    
214  Brick::~Brick()  Brick::~Brick()
215  {  {
     Paso_SystemMatrixPattern_free(m_pattern);  
     Paso_Connector_free(m_connector);  
     delete assembler;  
216  }  }
217    
218  string Brick::getDescription() const  string Brick::getDescription() const
# Line 412  void Brick::readNcGrid(escript::Data& ou Line 381  void Brick::readNcGrid(escript::Data& ou
381  #endif  #endif
382  }  }
383    
384    #ifdef USE_BOOSTIO
385    void Brick::readBinaryGridFromZipped(escript::Data& out, string filename,
386                               const ReaderParameters& params) const
387    {
388        // the mapping is not universally correct but should work on our
389        // supported platforms
390        switch (params.dataType) {
391            case DATATYPE_INT32:
392                readBinaryGridZippedImpl<int>(out, filename, params);
393                break;
394            case DATATYPE_FLOAT32:
395                readBinaryGridZippedImpl<float>(out, filename, params);
396                break;
397            case DATATYPE_FLOAT64:
398                readBinaryGridZippedImpl<double>(out, filename, params);
399                break;
400            default:
401                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
402        }
403    }
404    #endif
405    
406  void Brick::readBinaryGrid(escript::Data& out, string filename,  void Brick::readBinaryGrid(escript::Data& out, string filename,
407                             const ReaderParameters& params) const                             const ReaderParameters& params) const
408  {  {
# Line 548  void Brick::readBinaryGridImpl(escript:: Line 539  void Brick::readBinaryGridImpl(escript::
539      f.close();      f.close();
540  }  }
541    
542    #ifdef USE_BOOSTIO
543    template<typename ValueType>
544    void Brick::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
545                                   const ReaderParameters& params) const
546    {
547        // check destination function space
548        int myN0, myN1, myN2;
549        if (out.getFunctionSpace().getTypeCode() == Nodes) {
550            myN0 = m_NN[0];
551            myN1 = m_NN[1];
552            myN2 = m_NN[2];
553        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
554                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
555            myN0 = m_NE[0];
556            myN1 = m_NE[1];
557            myN2 = m_NE[2];
558        } else
559            throw RipleyException("readBinaryGridFromZipped(): invalid function space for output data object");
560    
561        if (params.first.size() != 3)
562            throw RipleyException("readBinaryGridFromZipped(): argument 'first' must have 3 entries");
563    
564        if (params.numValues.size() != 3)
565            throw RipleyException("readBinaryGridFromZipped(): argument 'numValues' must have 3 entries");
566    
567        if (params.multiplier.size() != 3)
568            throw RipleyException("readBinaryGridFromZipped(): argument 'multiplier' must have 3 entries");
569        for (size_t i=0; i<params.multiplier.size(); i++)
570            if (params.multiplier[i]<1)
571                throw RipleyException("readBinaryGridFromZipped(): all multipliers must be positive");
572    
573        // check file existence and size
574        ifstream f(filename.c_str(), ifstream::binary);
575        if (f.fail()) {
576            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
577        }
578        f.seekg(0, ios::end);
579        const int numComp = out.getDataPointSize();
580        int filesize = f.tellg();
581        f.seekg(0, ios::beg);
582        std::vector<char> compressed(filesize);
583        f.read((char*)&compressed[0], filesize);
584        f.close();
585        std::vector<char> decompressed = unzip(compressed);
586        filesize = decompressed.size();
587        const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
588        if (filesize < reqsize) {
589            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
590        }
591    
592        // check if this rank contributes anything
593        if (params.first[0] >= m_offset[0]+myN0 ||
594                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
595                params.first[1] >= m_offset[1]+myN1 ||
596                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
597                params.first[2] >= m_offset[2]+myN2 ||
598                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
599            return;
600        }
601    
602        // now determine how much this rank has to write
603    
604        // first coordinates in data object to write to
605        const int first0 = max(0, params.first[0]-m_offset[0]);
606        const int first1 = max(0, params.first[1]-m_offset[1]);
607        const int first2 = max(0, params.first[2]-m_offset[2]);
608        // indices to first value in file
609        const int idx0 = max(0, m_offset[0]-params.first[0]);
610        const int idx1 = max(0, m_offset[1]-params.first[1]);
611        const int idx2 = max(0, m_offset[2]-params.first[2]);
612        // number of values to read
613        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
614        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
615        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
616    
617        out.requireWrite();
618        vector<ValueType> values(num0*numComp);
619        const int dpp = out.getNumDataPointsPerSample();
620    
621        for (int z=0; z<num2; z++) {
622            for (int y=0; y<num1; y++) {
623                const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
624                                 +(idx2+z)*params.numValues[0]*params.numValues[1]);
625                memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
626                
627                for (int x=0; x<num0; x++) {
628                    const int baseIndex = first0+x*params.multiplier[0]
629                                         +(first1+y*params.multiplier[1])*myN0
630                                         +(first2+z*params.multiplier[2])*myN0*myN1;
631                    for (int m2=0; m2<params.multiplier[2]; m2++) {
632                        for (int m1=0; m1<params.multiplier[1]; m1++) {
633                            for (int m0=0; m0<params.multiplier[0]; m0++) {
634                                const int dataIndex = baseIndex+m0
635                                               +m1*myN0
636                                               +m2*myN0*myN1;
637                                double* dest = out.getSampleDataRW(dataIndex);
638                                for (int c=0; c<numComp; c++) {
639                                    ValueType val = values[x*numComp+c];
640    
641                                    if (params.byteOrder != BYTEORDER_NATIVE) {
642                                        char* cval = reinterpret_cast<char*>(&val);
643                                        // this will alter val!!
644                                        byte_swap32(cval);
645                                    }
646                                    if (!std::isnan(val)) {
647                                        for (int q=0; q<dpp; q++) {
648                                            *dest++ = static_cast<double>(val);
649                                        }
650                                    }
651                                }
652                            }
653                        }
654                    }
655                }
656            }
657        }
658    }
659    #endif
660    
661  void Brick::writeBinaryGrid(const escript::Data& in, string filename,  void Brick::writeBinaryGrid(const escript::Data& in, string filename,
662                              int byteOrder, int dataType) const                              int byteOrder, int dataType) const
663  {  {
# Line 1886  void Brick::assembleIntegrate(vector<dou Line 1996  void Brick::assembleIntegrate(vector<dou
1996  }  }
1997    
1998  //protected  //protected
1999  dim_t Brick::insertNeighbourNodes(IndexVector& index, index_t node) const  IndexVector Brick::getDiagonalIndices() const
2000  {  {
2001        IndexVector ret;
2002      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2003      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
     const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];  
     const int x=node%nDOF0;  
     const int y=node%(nDOF0*nDOF1)/nDOF0;  
     const int z=node/(nDOF0*nDOF1);  
     int num=0;  
     // loop through potential neighbours and add to index if positions are  
     // within bounds  
2004      for (int i2=-1; i2<2; i2++) {      for (int i2=-1; i2<2; i2++) {
2005          for (int i1=-1; i1<2; i1++) {          for (int i1=-1; i1<2; i1++) {
2006              for (int i0=-1; i0<2; i0++) {              for (int i0=-1; i0<2; i0++) {
2007                  // skip node itself                  ret.push_back(i2*nDOF0*nDOF1+i1*nDOF0+i0);
                 if (i0==0 && i1==0 && i2==0)  
                     continue;  
                 // location of neighbour node  
                 const int nx=x+i0;  
                 const int ny=y+i1;  
                 const int nz=z+i2;  
                 if (nx>=0 && ny>=0 && nz>=0  
                         && nx<nDOF0 && ny<nDOF1 && nz<nDOF2) {  
                     index.push_back(nz*nDOF0*nDOF1+ny*nDOF0+nx);  
                     num++;  
                 }  
2008              }              }
2009          }          }
2010      }      }
2011    
2012      return num;      return ret;
2013  }  }
2014    
2015  //protected  //protected
# Line 1947  void Brick::nodesToDOF(escript::Data& ou Line 2040  void Brick::nodesToDOF(escript::Data& ou
2040  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const
2041  {  {
2042      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
2043      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2044      // expand data object if necessary to be able to grab the whole data      // expand data object if necessary to be able to grab the whole data
2045      const_cast<escript::Data*>(&in)->expand();      const_cast<escript::Data*>(&in)->expand();
2046      Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));      coupler->startCollect(in.getDataRO());
2047    
2048      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2049      out.requireWrite();      out.requireWrite();
2050      const double* buffer = Paso_Coupler_finishCollect(coupler);      const double* buffer = coupler->finishCollect();
2051    
2052  #pragma omp parallel for  #pragma omp parallel for
2053      for (index_t i=0; i<getNumNodes(); i++) {      for (index_t i=0; i<getNumNodes(); i++) {
# Line 1963  void Brick::dofToNodes(escript::Data& ou Line 2056  void Brick::dofToNodes(escript::Data& ou
2056                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
2057          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
2058      }      }
     Paso_Coupler_free(coupler);  
2059  }  }
2060    
2061  //private  //private
# Line 1984  void Brick::populateSampleIds() Line 2076  void Brick::populateSampleIds()
2076          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2077      }      }
2078      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2079      m_nodeId.resize(getNumNodes());      
2080      m_dofId.resize(numDOF);      try {
2081      m_elementId.resize(getNumElements());          m_nodeId.resize(getNumNodes());
2082            m_dofId.resize(numDOF);
2083            m_elementId.resize(getNumElements());
2084        } catch (const std::length_error& le) {
2085            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2086        }
2087        
2088      // populate face element counts      // populate face element counts
2089      //left      //left
2090      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2231  void Brick::createPattern() Line 2328  void Brick::createPattern()
2328      RankVector neighbour;      RankVector neighbour;
2329      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
2330      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
2331      int numShared=0;      int numShared=0, expectedShared=0;;
2332      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
2333      const int y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];      const int y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
2334      const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
# Line 2245  void Brick::createPattern() Line 2342  void Brick::createPattern()
2342                  const int nx=x+i0;                  const int nx=x+i0;
2343                  const int ny=y+i1;                  const int ny=y+i1;
2344                  const int nz=z+i2;                  const int nz=z+i2;
2345                    if (!(nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2])) {
2346                        continue;
2347                    }
2348                    if (i0==0 && i1==0)
2349                        expectedShared += nDOF0*nDOF1;
2350                    else if (i0==0 && i2==0)
2351                        expectedShared += nDOF0*nDOF2;
2352                    else if (i1==0 && i2==0)
2353                        expectedShared += nDOF1*nDOF2;
2354                    else if (i0==0)
2355                        expectedShared += nDOF0;
2356                    else if (i1==0)
2357                        expectedShared += nDOF1;
2358                    else if (i2==0)
2359                        expectedShared += nDOF2;
2360                    else
2361                        expectedShared++;
2362                }
2363            }
2364        }
2365        
2366        vector<IndexVector> rowIndices(expectedShared);
2367        
2368        for (int i2=-1; i2<2; i2++) {
2369            for (int i1=-1; i1<2; i1++) {
2370                for (int i0=-1; i0<2; i0++) {
2371                    // skip this rank
2372                    if (i0==0 && i1==0 && i2==0)
2373                        continue;
2374                    // location of neighbour rank
2375                    const int nx=x+i0;
2376                    const int ny=y+i1;
2377                    const int nz=z+i2;
2378                  if (nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2]) {                  if (nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2]) {
2379                      neighbour.push_back(nz*m_NX[0]*m_NX[1]+ny*m_NX[0]+nx);                      neighbour.push_back(nz*m_NX[0]*m_NX[1]+ny*m_NX[0]+nx);
2380                      if (i0==0 && i1==0) {                      if (i0==0 && i1==0) {
# Line 2260  void Brick::createPattern() Line 2390  void Brick::createPattern()
2390                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2391                                  if (j>0) {                                  if (j>0) {
2392                                      if (i>0)                                      if (i>0)
2393                                          colIndices[firstDOF+j-1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0, numShared);
2394                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2395                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2396                                          colIndices[firstDOF+j-1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0, numShared);
2397                                  }                                  }
2398                                  if (i>0)                                  if (i>0)
2399                                      colIndices[firstDOF+j-nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0, numShared);
2400                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2401                                  if (i<nDOF1-1)                                  if (i<nDOF1-1)
2402                                      colIndices[firstDOF+j+nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0, numShared);
2403                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2404                                      if (i>0)                                      if (i>0)
2405                                          colIndices[firstDOF+j+1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0, numShared);
2406                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2407                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2408                                          colIndices[firstDOF+j+1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0, numShared);
2409                                  }                                  }
2410                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2411                              }                              }
# Line 2294  void Brick::createPattern() Line 2424  void Brick::createPattern()
2424                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2425                                  if (j>0) {                                  if (j>0) {
2426                                      if (i>0)                                      if (i>0)
2427                                          colIndices[firstDOF+j-1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0*nDOF1, numShared);
2428                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2429                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2430                                          colIndices[firstDOF+j-1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0*nDOF1, numShared);
2431                                  }                                  }
2432                                  if (i>0)                                  if (i>0)
2433                                      colIndices[firstDOF+j-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0*nDOF1, numShared);
2434                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2435                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2436                                      colIndices[firstDOF+j+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0*nDOF1, numShared);
2437                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2438                                      if (i>0)                                      if (i>0)
2439                                          colIndices[firstDOF+j+1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0*nDOF1, numShared);
2440                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2441                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2442                                          colIndices[firstDOF+j+1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0*nDOF1, numShared);
2443                                  }                                  }
2444                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2445                              }                              }
# Line 2328  void Brick::createPattern() Line 2458  void Brick::createPattern()
2458                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2459                                  if (j>0) {                                  if (j>0) {
2460                                      if (i>0)                                      if (i>0)
2461                                          colIndices[firstDOF+(j-1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0-nDOF0*nDOF1, numShared);
2462                                      colIndices[firstDOF+(j-1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0, numShared);
2463                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2464                                          colIndices[firstDOF+(j-1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0+nDOF0*nDOF1, numShared);
2465                                  }                                  }
2466                                  if (i>0)                                  if (i>0)
2467                                      colIndices[firstDOF+j*nDOF0-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0-nDOF0*nDOF1, numShared);
2468                                  colIndices[firstDOF+j*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0, numShared);
2469                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2470                                      colIndices[firstDOF+j*nDOF0+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0+nDOF0*nDOF1, numShared);
2471                                  if (j<nDOF1-1) {                                  if (j<nDOF1-1) {
2472                                      if (i>0)                                      if (i>0)
2473                                          colIndices[firstDOF+(j+1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0-nDOF0*nDOF1, numShared);
2474                                      colIndices[firstDOF+(j+1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0, numShared);
2475                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2476                                          colIndices[firstDOF+(j+1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0+nDOF0*nDOF1, numShared);
2477                                  }                                  }
2478                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2479                              }                              }
# Line 2359  void Brick::createPattern() Line 2489  void Brick::createPattern()
2489                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
2490                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2491                              if (i>0)                              if (i>0)
2492                                  colIndices[firstDOF+i-1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
2493                              colIndices[firstDOF+i].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
2494                              if (i<nDOF0-1)                              if (i<nDOF0-1)
2495                                  colIndices[firstDOF+i+1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
2496                              m_dofMap[firstNode+i]=numDOF+numShared;                              m_dofMap[firstNode+i]=numDOF+numShared;
2497                          }                          }
2498                      } else if (i1==0) {                      } else if (i1==0) {
# Line 2377  void Brick::createPattern() Line 2507  void Brick::createPattern()
2507                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
2508                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2509                              if (i>0)                              if (i>0)
2510                                  colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
2511                              colIndices[firstDOF+i*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
2512                              if (i<nDOF1-1)                              if (i<nDOF1-1)
2513                                  colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
2514                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2515                          }                          }
2516                      } else if (i2==0) {                      } else if (i2==0) {
# Line 2395  void Brick::createPattern() Line 2525  void Brick::createPattern()
2525                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2526                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2527                              if (i>0)                              if (i>0)
2528                                  colIndices[firstDOF+(i-1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0*nDOF1, numShared);
2529                              colIndices[firstDOF+i*nDOF0*nDOF1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0*nDOF1, numShared);
2530                              if (i<nDOF2-1)                              if (i<nDOF2-1)
2531                                  colIndices[firstDOF+(i+1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0*nDOF1, numShared);
2532                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2533                          }                          }
2534                      } else {                      } else {
# Line 2412  void Brick::createPattern() Line 2542  void Brick::createPattern()
2542                          offsetInShared.push_back(offsetInShared.back()+1);                          offsetInShared.push_back(offsetInShared.back()+1);
2543                          sendShared.push_back(dof);                          sendShared.push_back(dof);
2544                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
2545                          colIndices[dof].push_back(numShared);                          doublyLink(colIndices, rowIndices, dof, numShared);
2546                          m_dofMap[node]=numDOF+numShared;                          m_dofMap[node]=numDOF+numShared;
2547                          ++numShared;                          ++numShared;
2548                      }                      }
# Line 2421  void Brick::createPattern() Line 2551  void Brick::createPattern()
2551          }          }
2552      }      }
2553    
2554    #pragma omp parallel for
2555        for (int i = 0; i < numShared; i++) {
2556            std::sort(rowIndices[i].begin(), rowIndices[i].end());
2557        }
2558    
2559      // create connector      // create connector
2560      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
2561              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2562              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2563      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2564              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2565              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2566      m_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
     Paso_SharedComponents_free(snd_shcomp);  
     Paso_SharedComponents_free(rcv_shcomp);  
   
     // create main and couple blocks  
     Paso_Pattern *mainPattern = createMainPattern();  
     Paso_Pattern *colPattern, *rowPattern;  
     createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);  
   
     // allocate paso distribution  
     Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,  
             const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);  
   
     // finally create the system matrix  
     m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,  
             distribution, distribution, mainPattern, colPattern, rowPattern,  
             m_connector, m_connector);  
   
     Paso_Distribution_free(distribution);  
2567    
2568      // useful debug output      // useful debug output
2569      /*      /*
# Line 2505  void Brick::createPattern() Line 2622  void Brick::createPattern()
2622          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2623      }      }
2624      */      */
   
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colPattern);  
     Paso_Pattern_free(rowPattern);  
2625  }  }
2626    
2627  //private  //private
2628  void Brick::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,  void Brick::addToMatrixAndRHS(SystemMatrix* S, escript::Data& F,
2629           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2630           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2631  {  {
# Line 2536  void Brick::addToMatrixAndRHS(Paso_Syste Line 2649  void Brick::addToMatrixAndRHS(Paso_Syste
2649          }          }
2650      }      }
2651      if (addS) {      if (addS) {
2652          addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);          S->add(rowIndex, EM_S);
2653      }      }
2654  }  }
2655    
# Line 2861  void Brick::interpolateNodesOnFaces(escr Line 2974  void Brick::interpolateNodesOnFaces(escr
2974    
2975  namespace  namespace
2976  {  {
2977      // Calculates a guassian blur colvolution matrix for 3D      // Calculates a gaussian blur convolution matrix for 3D
2978      // See wiki article on the subject      // See wiki article on the subject
2979      double* get3DGauss(unsigned radius, double sigma)      double* get3DGauss(unsigned radius, double sigma)
2980      {      {
2981          double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];          double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];
2982          double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);          double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);
2983      double total=0;          double total=0;
2984      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
2985      for (int z=-r;z<=r;++z)          for (int z=-r;z<=r;++z)
2986      {          {
2987          for (int y=-r;y<=r;++y)              for (int y=-r;y<=r;++y)
2988          {              {
2989          for (int x=-r;x<=r;++x)                  for (int x=-r;x<=r;++x)
2990          {                          {        
2991              arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));                      arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));
2992              total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];                      total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];    
2993          }                  }
2994          }              }
2995      }          }
2996      double invtotal=1/total;          double invtotal=1/total;
2997      for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)          for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
2998      {          {
2999          arr[p]*=invtotal;              arr[p]*=invtotal;
3000      }          }
3001      return arr;          return arr;
3002      }      }
3003            
3004      // applies conv to source to get a point.      // applies conv to source to get a point.
# Line 2893  namespace Line 3006  namespace
3006      double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)      double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)
3007      {      {
3008          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3009      size_t sbase=bx+by*width+bz*width*height;          size_t sbase=bx+by*width+bz*width*height;
3010      double result=0;          double result=0;
3011      for (int z=0;z<2*radius+1;++z)          for (int z=0;z<2*radius+1;++z)
3012      {          {
3013          for (int y=0;y<2*radius+1;++y)              for (int y=0;y<2*radius+1;++y)
3014          {                  {    
3015          for (int x=0;x<2*radius+1;++x)                  for (int x=0;x<2*radius+1;++x)
3016          {                  {
3017              result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];                      result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];
3018          }                  }
3019          }              }
3020      }          }
3021      // use this line for "pass-though" (return the centre point value)          // use this line for "pass-though" (return the centre point value)
3022  //  return source[sbase+(radius)+(radius)*width+(radius)*width*height];  //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3023      return result;                return result;      
3024      }      }
3025  }  }
3026    
# Line 2921  escript::Data Brick::randomFill(const es Line 3034  escript::Data Brick::randomFill(const es
3034      int numvals=escript::DataTypes::noValues(shape);      int numvals=escript::DataTypes::noValues(shape);
3035      if (len(filter)>0 && (numvals!=1))      if (len(filter)>0 && (numvals!=1))
3036      {      {
3037      throw RipleyException("Ripley only supports filters for scalar data.");          throw RipleyException("Ripley only supports filters for scalar data.");
3038      }      }
3039      escript::Data res=randomFillWorker(shape, seed, filter);      escript::Data res=randomFillWorker(shape, seed, filter);
3040      if (res.getFunctionSpace()!=what)      if (res.getFunctionSpace()!=what)
3041      {      {
3042      escript::Data r=escript::Data(res, what);          escript::Data r=escript::Data(res, what);
3043      return r;          return r;
3044      }      }
3045      return res;      return res;
3046  }  }
# Line 2971  that ripley has. Line 3084  that ripley has.
3084  */  */
3085  escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const  escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
3086  {  {
3087      if (m_numDim!=3)      unsigned int radius=0;  // these are only used by gaussian
     {  
         throw RipleyException("Brick must be 3D.");  
     }  
       
     unsigned int radius=0;  // these are only used by gaussian  
3088      double sigma=0.5;      double sigma=0.5;
3089            
3090      unsigned int numvals=escript::DataTypes::noValues(shape);      unsigned int numvals=escript::DataTypes::noValues(shape);
3091            
3092      if (len(filter)==0)      if (len(filter)==0)
3093      {      {
3094      // nothing special required here yet      // nothing special required here yet
3095      }      }
3096      else if (len(filter)==3)      else if (len(filter)==3)
3097      {      {
3098      boost::python::extract<string> ex(filter[0]);          boost::python::extract<string> ex(filter[0]);
3099      if (!ex.check() || (ex()!="gaussian"))          if (!ex.check() || (ex()!="gaussian"))
3100      {          {
3101          throw RipleyException("Unsupported random filter for Brick.");              throw RipleyException("Unsupported random filter for Brick.");
3102      }          }
3103      boost::python::extract<unsigned int> ex1(filter[1]);          boost::python::extract<unsigned int> ex1(filter[1]);
3104      if (!ex1.check())          if (!ex1.check())
3105      {          {
3106          throw RipleyException("Radius of gaussian filter must be a positive integer.");              throw RipleyException("Radius of gaussian filter must be a positive integer.");
3107      }          }
3108      radius=ex1();          radius=ex1();
3109      sigma=0.5;          sigma=0.5;
3110      boost::python::extract<double> ex2(filter[2]);          boost::python::extract<double> ex2(filter[2]);
3111      if (!ex2.check() || (sigma=ex2())<=0)          if (!ex2.check() || (sigma=ex2())<=0)
3112      {          {
3113          throw RipleyException("Sigma must be a postive floating point number.");              throw RipleyException("Sigma must be a positive floating point number.");
3114      }                      }            
3115      }      }
3116      else      else
3117      {      {
3118          throw RipleyException("Unsupported random filter");          throw RipleyException("Unsupported random filter");
3119      }      }
       
       
3120    
3121            // number of points in the internal region
3122      size_t internal[3];      // that is, the ones we need smoothed versions of
3123      internal[0]=m_NE[0]+1;  // number of points in the internal region      const dim_t internal[3] = { m_NN[0], m_NN[1], m_NN[2] };
     internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of  
     internal[2]=m_NE[2]+1;  // that is, the ones we need smoothed versions of  
3124      size_t ext[3];      size_t ext[3];
3125      ext[0]=internal[0]+2*radius;    // includes points we need as input      ext[0]=(size_t)internal[0]+2*radius;  // includes points we need as input
3126      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3127      ext[2]=internal[2]+2*radius;    // for smoothing      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3128            
3129      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3130      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
# Line 3045  escript::Data Brick::randomFillWorker(co Line 3149  escript::Data Brick::randomFillWorker(co
3149        
3150      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3151      {      {
3152      // since the dimensions are equal for all ranks, this exception      // since the dimensions are equal for all ranks, this exception
3153      // will be thrown on all ranks      // will be thrown on all ranks
3154      throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");      throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3155    
3156      }      }
3157      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
# Line 3063  escript::Data Brick::randomFillWorker(co Line 3167  escript::Data Brick::randomFillWorker(co
3167  #ifdef ESYS_MPI      #ifdef ESYS_MPI    
3168      basex=X*m_gNE[0]/m_NX[0];      basex=X*m_gNE[0]/m_NX[0];
3169      basey=Y*m_gNE[1]/m_NX[1];      basey=Y*m_gNE[1]/m_NX[1];
3170      basez=Z*m_gNE[2]/m_NX[2];          basez=Z*m_gNE[2]/m_NX[2];
 #endif      
     if (seed==0)  
     {  
     seed=2; // since we are using the seed parameter as the spacing and 0 spacing causes an exception  
     }  
     esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, seed, basex, basey, basez);  
 */  
       
3171            
3172  /*  cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
 cout << "Pattern:\n";      
 for (int i=0;i<ext[0]*ext[1]*ext[2];)  
 {  
     cout << src[i++] << " ";  
     if (i%ext[0]==0)  
       cout << "\n";  
     if (i%(ext[0]*ext[1])==0)  
       cout << "\n";  
 }*/  
3173            
3174      #endif    
3175        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3176    */
3177    
3178  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3179    
3180    
3181    
3182      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);
3183      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3184          // there is an overlap of two points both of which need to have "radius" points on either side.          // there is an overlap of two points both of which need to have "radius" points on either side.
3185            
3186      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets
3187      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;  
3188      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3189            
3190      Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);          Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);    
3191            
3192      MPI_Request reqs[50];       // a non-tight upper bound on how many we need      MPI_Request reqs[50];       // a non-tight upper bound on how many we need
3193      MPI_Status stats[50];      MPI_Status stats[50];
3194      short rused=0;      short rused=0;
3195            
# Line 3112  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3202  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3202            
3203      block.copyAllToBuffer(src);      block.copyAllToBuffer(src);
3204            
       
3205      int comserr=0;          int comserr=0;    
3206      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
3207      {      {
3208      message& m=incoms[i];          message& m=incoms[i];
3209      comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3210      block.setUsed(m.destbuffid);          block.setUsed(m.destbuffid);
3211      }      }
3212    
3213      for (size_t i=0;i<outcoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3214      {      {
3215      message& m=outcoms[i];          message& m=outcoms[i];
3216      comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3217      }          }    
3218            
3219      if (!comserr)      if (!comserr)
# Line 3134  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3223  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3223    
3224      if (comserr)      if (comserr)
3225      {      {
3226      // Yes this is throwing an exception as a result of an MPI error.      // Yes this is throwing an exception as a result of an MPI error.
3227      // and no we don't inform the other ranks that we are doing this.      // and no we don't inform the other ranks that we are doing this.
3228      // however, we have no reason to believe coms work at this point anyway      // however, we have no reason to believe coms work at this point anyway
3229          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
3230      }      }
3231            
3232      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3233            
       
       
   
3234  #endif      #endif    
3235            
3236  /*          if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
 cout << "Pattern (post transfer):\n";      
 for (int i=0;i<ext[0]*ext[1]*ext[2];)  
 {  
     cout << src[i++] << " ";  
     if (i%ext[0]==0)  
       cout << "\n";  
     if (i%(ext[0]*ext[1])==0)  
       cout << "\n";  
 }   */  
       
     if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe  
3237      {      {
3238                
3239      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3240      escript::Data resdat(0, shape, fs , true);          escript::Data resdat(0, shape, fs , true);
3241      // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3242      escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3243            
3244            
3245      // now we need to copy values over          // now we need to copy values over
3246      for (size_t z=0;z<(internal[2]);++z)          for (size_t z=0;z<(internal[2]);++z)
3247      {          {
3248          for (size_t y=0;y<(internal[1]);++y)                  for (size_t y=0;y<(internal[1]);++y)    
3249          {              {
3250          for (size_t x=0;x<(internal[0]);++x)                  for (size_t x=0;x<(internal[0]);++x)
3251          {                  {
3252              for (unsigned int i=0;i<numvals;++i)                      for (unsigned int i=0;i<numvals;++i)
3253              {                      {
3254                  dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];                          dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3255              }                      }
3256          }                  }
3257          }              }
3258      }            }  
3259      delete[] src;          delete[] src;
3260      return resdat;                return resdat;      
3261      }      }
3262      else        // filter enabled        else        // filter enabled  
3263      {      {
3264            
3265      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3266      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3267      // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3268      escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3269      double* convolution=get3DGauss(radius, sigma);          double* convolution=get3DGauss(radius, sigma);
3270    
3271  // cout << "Convolution matrix\n";          for (size_t z=0;z<(internal[2]);++z)
3272  // size_t di=(radius*2+1);          {
3273  // double ctot=0;              for (size_t y=0;y<(internal[1]);++y)    
3274  // for (int i=0;i<di*di*di;++i)              {
3275  // {                  for (size_t x=0;x<(internal[0]);++x)
3276  //     cout << convolution[i] << " ";                  {    
3277  //     ctot+=convolution[i];                      dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3278  //     if ((i+1)%di==0)              
3279  //     {                  }
3280  //  cout << "\n";              }
3281  //     }          }
3282  //     if ((i+1)%(di*di)==0)      
3283  //     {          delete[] convolution;
3284  //  cout << "\n";          delete[] src;
3285  //     }          return resdat;
 // }  
 //  
 // cout << "Sum of matrix is = " << ctot << endl;  
   
     for (size_t z=0;z<(internal[2]);++z)  
     {  
         for (size_t y=0;y<(internal[1]);++y)      
         {  
         for (size_t x=0;x<(internal[0]);++x)  
         {      
             dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);  
               
         }  
         }  
     }  
       
 // cout << "\nResulting matrix:\n";  
 //     for (size_t z=0;z<(internal[2]);++z)  
 //     {  
 //  for (size_t y=0;y<(internal[1]);++y)      
 //  {  
 //      for (size_t x=0;x<(internal[0]);++x)  
 //      {      
 //      cout << dv[x+y*(internal[0])+z*internal[0]*internal[1]] << " ";  
 //      }  
 //      cout << endl;  
 //  }  
 //  cout << endl;  
 //     }  
   
       
     delete[] convolution;  
     delete[] src;  
     return resdat;  
3286            
3287      }      }
3288  }  }
3289    
3290    int Brick::findNode(const double *coords) const
3291    {
   
   
   
 int Brick::findNode(const double *coords) const {  
3292      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3293      //is the found element even owned by this rank      //is the found element even owned by this rank
3294      // (inside owned or shared elements but will map to an owned element)      // (inside owned or shared elements but will map to an owned element)
# Line 3268  int Brick::findNode(const double *coords Line 3305  int Brick::findNode(const double *coords
3305      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3306      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3307      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3308        
3309        //check if the point is even inside the domain
3310        if (x < 0 || y < 0 || z < 0
3311                || x > m_length[0] || y > m_length[1] || z > m_length[2])
3312            return NOT_MINE;
3313            
3314      // distance in elements      // distance in elements
3315      int ex = (int) floor(x / m_dx[0]);      int ex = (int) floor(x / m_dx[0]);
3316      int ey = (int) floor(y / m_dx[1]);      int ey = (int) floor(y / m_dx[1]);
# Line 3300  int Brick::findNode(const double *coords Line 3343  int Brick::findNode(const double *coords
3343      return closest;      return closest;
3344  }  }
3345    
3346  void Brick::setAssembler(std::string type, std::map<std::string,  Assembler_ptr Brick::createAssembler(std::string type, std::map<std::string,
3347          escript::Data> constants) {          escript::Data> constants) const
3348      if (type.compare("WaveAssembler") == 0) {  {
3349          delete assembler;      if (type.compare("DefaultAssembler") == 0) {
3350          assembler = new WaveAssembler3D(this, m_dx, m_NX, m_NE, m_NN, constants);          return Assembler_ptr(new DefaultAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
3351        } else if (type.compare("WaveAssembler") == 0) {
3352            return Assembler_ptr(new WaveAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN, constants));
3353        } else if (type.compare("LameAssembler") == 0) {
3354            return Assembler_ptr(new LameAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
3355      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3356          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Brick does not support the"
3357                                  " requested assembler");                                  " requested assembler");
3358      }      }
3359  }  }

Legend:
Removed from v.4687  
changed lines
  Added in v.4949

  ViewVC Help
Powered by ViewVC 1.1.26