/[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 4653 by jfenwick, Wed Feb 5 05:20:41 2014 UTC branches/diaplayground/ripley/src/Brick.cpp revision 4988 by caltinay, Wed Jun 4 01:15:10 2014 UTC
# Line 1  Line 1 
1    
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2013 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
# Line 9  Line 9 
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  * Development since 2012 by School of Earth Sciences  * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
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 42  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 55  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 229  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 240  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 411  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 460  void Brick::readBinaryGridImpl(escript:: Line 452  void Brick::readBinaryGridImpl(escript::
452      for (size_t i=0; i<params.multiplier.size(); i++)      for (size_t i=0; i<params.multiplier.size(); i++)
453          if (params.multiplier[i]<1)          if (params.multiplier[i]<1)
454              throw RipleyException("readBinaryGrid(): all multipliers must be positive");              throw RipleyException("readBinaryGrid(): all multipliers must be positive");
455        if (params.reverse[0] != 0 || params.reverse[1] != 0)
456            throw RipleyException("readBinaryGrid(): reversing only supported in Z-direction currently");
457    
458      // check file existence and size      // check file existence and size
459      ifstream f(filename.c_str(), ifstream::binary);      ifstream f(filename.c_str(), ifstream::binary);
# Line 492  void Brick::readBinaryGridImpl(escript:: Line 486  void Brick::readBinaryGridImpl(escript::
486      const int first0 = max(0, params.first[0]-m_offset[0]);      const int first0 = max(0, params.first[0]-m_offset[0]);
487      const int first1 = max(0, params.first[1]-m_offset[1]);      const int first1 = max(0, params.first[1]-m_offset[1]);
488      const int first2 = max(0, params.first[2]-m_offset[2]);      const int first2 = max(0, params.first[2]-m_offset[2]);
489        // indices to first value in file (not accounting for reverse yet)
490        int idx0 = max(0, m_offset[0]-params.first[0]);
491        int idx1 = max(0, m_offset[1]-params.first[1]);
492        int idx2 = max(0, m_offset[2]-params.first[2]);
493        // number of values to read
494        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
495        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
496        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
497    
498        // make sure we read the right block if going backwards through file
499        if (params.reverse[2])
500            idx2 = params.numValues[2]-idx2-1;
501    
502        // helpers for reversing
503        const int z_mult = (params.reverse[2] ? -1 : 1);
504    
505        out.requireWrite();
506        vector<ValueType> values(num0*numComp);
507        const int dpp = out.getNumDataPointsPerSample();
508    
509        for (int z=0; z<num2; z++) {
510            for (int y=0; y<num1; y++) {
511                const int fileofs = numComp*(idx0 +
512                                    (idx1+y)*params.numValues[0] +
513                                    (idx2+z_mult*z)*params.numValues[0]*params.numValues[1]);
514                f.seekg(fileofs*sizeof(ValueType));
515                f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
516    
517                for (int x=0; x<num0; x++) {
518                    const int baseIndex = first0+x*params.multiplier[0]
519                                         +(first1+y*params.multiplier[1])*myN0
520                                         +(first2+z*params.multiplier[2])*myN0*myN1;
521                    for (int m2=0; m2<params.multiplier[2]; m2++) {
522                        for (int m1=0; m1<params.multiplier[1]; m1++) {
523                            for (int m0=0; m0<params.multiplier[0]; m0++) {
524                                const int dataIndex = baseIndex+m0
525                                               +m1*myN0
526                                               +m2*myN0*myN1;
527                                double* dest = out.getSampleDataRW(dataIndex);
528                                for (int c=0; c<numComp; c++) {
529                                    ValueType val = values[x*numComp+c];
530    
531                                    if (params.byteOrder != BYTEORDER_NATIVE) {
532                                        char* cval = reinterpret_cast<char*>(&val);
533                                        // this will alter val!!
534                                        byte_swap32(cval);
535                                    }
536                                    if (!std::isnan(val)) {
537                                        for (int q=0; q<dpp; q++) {
538                                            *dest++ = static_cast<double>(val);
539                                        }
540                                    }
541                                }
542                            }
543                        }
544                    }
545                }
546            }
547        }
548    
549        f.close();
550    }
551    
552    #ifdef USE_BOOSTIO
553    template<typename ValueType>
554    void Brick::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
555                                   const ReaderParameters& params) const
556    {
557        // check destination function space
558        int myN0, myN1, myN2;
559        if (out.getFunctionSpace().getTypeCode() == Nodes) {
560            myN0 = m_NN[0];
561            myN1 = m_NN[1];
562            myN2 = m_NN[2];
563        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
564                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
565            myN0 = m_NE[0];
566            myN1 = m_NE[1];
567            myN2 = m_NE[2];
568        } else
569            throw RipleyException("readBinaryGridFromZipped(): invalid function space for output data object");
570    
571        if (params.first.size() != 3)
572            throw RipleyException("readBinaryGridFromZipped(): argument 'first' must have 3 entries");
573    
574        if (params.numValues.size() != 3)
575            throw RipleyException("readBinaryGridFromZipped(): argument 'numValues' must have 3 entries");
576    
577        if (params.multiplier.size() != 3)
578            throw RipleyException("readBinaryGridFromZipped(): argument 'multiplier' must have 3 entries");
579        for (size_t i=0; i<params.multiplier.size(); i++)
580            if (params.multiplier[i]<1)
581                throw RipleyException("readBinaryGridFromZipped(): all multipliers must be positive");
582    
583        // check file existence and size
584        ifstream f(filename.c_str(), ifstream::binary);
585        if (f.fail()) {
586            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
587        }
588        f.seekg(0, ios::end);
589        const int numComp = out.getDataPointSize();
590        int filesize = f.tellg();
591        f.seekg(0, ios::beg);
592        std::vector<char> compressed(filesize);
593        f.read((char*)&compressed[0], filesize);
594        f.close();
595        std::vector<char> decompressed = unzip(compressed);
596        filesize = decompressed.size();
597        const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
598        if (filesize < reqsize) {
599            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
600        }
601    
602        // check if this rank contributes anything
603        if (params.first[0] >= m_offset[0]+myN0 ||
604                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
605                params.first[1] >= m_offset[1]+myN1 ||
606                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
607                params.first[2] >= m_offset[2]+myN2 ||
608                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
609            return;
610        }
611    
612        // now determine how much this rank has to write
613    
614        // first coordinates in data object to write to
615        const int first0 = max(0, params.first[0]-m_offset[0]);
616        const int first1 = max(0, params.first[1]-m_offset[1]);
617        const int first2 = max(0, params.first[2]-m_offset[2]);
618      // indices to first value in file      // indices to first value in file
619      const int idx0 = max(0, m_offset[0]-params.first[0]);      const int idx0 = max(0, m_offset[0]-params.first[0]);
620      const int idx1 = max(0, m_offset[1]-params.first[1]);      const int idx1 = max(0, m_offset[1]-params.first[1]);
# Line 509  void Brick::readBinaryGridImpl(escript:: Line 632  void Brick::readBinaryGridImpl(escript::
632          for (int y=0; y<num1; y++) {          for (int y=0; y<num1; y++) {
633              const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]              const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
634                               +(idx2+z)*params.numValues[0]*params.numValues[1]);                               +(idx2+z)*params.numValues[0]*params.numValues[1]);
635              f.seekg(fileofs*sizeof(ValueType));              memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
636              f.read((char*)&values[0], num0*numComp*sizeof(ValueType));              
   
637              for (int x=0; x<num0; x++) {              for (int x=0; x<num0; x++) {
638                  const int baseIndex = first0+x*params.multiplier[0]                  const int baseIndex = first0+x*params.multiplier[0]
639                                       +(first1+y*params.multiplier[1])*myN0                                       +(first1+y*params.multiplier[1])*myN0
# Line 543  void Brick::readBinaryGridImpl(escript:: Line 665  void Brick::readBinaryGridImpl(escript::
665              }              }
666          }          }
667      }      }
   
     f.close();  
668  }  }
669    #endif
670    
671  void Brick::writeBinaryGrid(const escript::Data& in, string filename,  void Brick::writeBinaryGrid(const escript::Data& in, string filename,
672                              int byteOrder, int dataType) const                              int byteOrder, int dataType) const
# Line 786  const int* Brick::borrowSampleReferenceI Line 907  const int* Brick::borrowSampleReferenceI
907          case FaceElements:          case FaceElements:
908          case ReducedFaceElements:          case ReducedFaceElements:
909              return &m_faceId[0];              return &m_faceId[0];
910            case Points:
911                return &m_diracPointNodeIDs[0];
912          default:          default:
913              break;              break;
914      }      }
# Line 1883  void Brick::assembleIntegrate(vector<dou Line 2006  void Brick::assembleIntegrate(vector<dou
2006  }  }
2007    
2008  //protected  //protected
2009  dim_t Brick::insertNeighbourNodes(IndexVector& index, index_t node) const  IndexVector Brick::getDiagonalIndices() const
2010  {  {
2011        IndexVector ret;
2012      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2013      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  
2014      for (int i2=-1; i2<2; i2++) {      for (int i2=-1; i2<2; i2++) {
2015          for (int i1=-1; i1<2; i1++) {          for (int i1=-1; i1<2; i1++) {
2016              for (int i0=-1; i0<2; i0++) {              for (int i0=-1; i0<2; i0++) {
2017                  // 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++;  
                 }  
2018              }              }
2019          }          }
2020      }      }
2021    
2022      return num;      return ret;
2023  }  }
2024    
2025  //protected  //protected
# Line 1944  void Brick::nodesToDOF(escript::Data& ou Line 2050  void Brick::nodesToDOF(escript::Data& ou
2050  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const
2051  {  {
2052      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
2053      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2054      // 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
2055      const_cast<escript::Data*>(&in)->expand();      const_cast<escript::Data*>(&in)->expand();
2056      Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));      coupler->startCollect(in.getDataRO());
2057    
2058      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2059      out.requireWrite();      out.requireWrite();
2060      const double* buffer = Paso_Coupler_finishCollect(coupler);      const double* buffer = coupler->finishCollect();
2061    
2062  #pragma omp parallel for  #pragma omp parallel for
2063      for (index_t i=0; i<getNumNodes(); i++) {      for (index_t i=0; i<getNumNodes(); i++) {
# Line 1960  void Brick::dofToNodes(escript::Data& ou Line 2066  void Brick::dofToNodes(escript::Data& ou
2066                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
2067          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
2068      }      }
     Paso_Coupler_free(coupler);  
2069  }  }
2070    
2071  //private  //private
# Line 1981  void Brick::populateSampleIds() Line 2086  void Brick::populateSampleIds()
2086          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2087      }      }
2088      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2089      m_nodeId.resize(getNumNodes());      
2090      m_dofId.resize(numDOF);      try {
2091      m_elementId.resize(getNumElements());          m_nodeId.resize(getNumNodes());
2092            m_dofId.resize(numDOF);
2093            m_elementId.resize(getNumElements());
2094        } catch (const std::length_error& le) {
2095            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2096        }
2097        
2098      // populate face element counts      // populate face element counts
2099      //left      //left
2100      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2228  void Brick::createPattern() Line 2338  void Brick::createPattern()
2338      RankVector neighbour;      RankVector neighbour;
2339      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
2340      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
2341      int numShared=0;      int numShared=0, expectedShared=0;;
2342      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
2343      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];
2344      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 2242  void Brick::createPattern() Line 2352  void Brick::createPattern()
2352                  const int nx=x+i0;                  const int nx=x+i0;
2353                  const int ny=y+i1;                  const int ny=y+i1;
2354                  const int nz=z+i2;                  const int nz=z+i2;
2355                    if (!(nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2])) {
2356                        continue;
2357                    }
2358                    if (i0==0 && i1==0)
2359                        expectedShared += nDOF0*nDOF1;
2360                    else if (i0==0 && i2==0)
2361                        expectedShared += nDOF0*nDOF2;
2362                    else if (i1==0 && i2==0)
2363                        expectedShared += nDOF1*nDOF2;
2364                    else if (i0==0)
2365                        expectedShared += nDOF0;
2366                    else if (i1==0)
2367                        expectedShared += nDOF1;
2368                    else if (i2==0)
2369                        expectedShared += nDOF2;
2370                    else
2371                        expectedShared++;
2372                }
2373            }
2374        }
2375        
2376        vector<IndexVector> rowIndices(expectedShared);
2377        
2378        for (int i2=-1; i2<2; i2++) {
2379            for (int i1=-1; i1<2; i1++) {
2380                for (int i0=-1; i0<2; i0++) {
2381                    // skip this rank
2382                    if (i0==0 && i1==0 && i2==0)
2383                        continue;
2384                    // location of neighbour rank
2385                    const int nx=x+i0;
2386                    const int ny=y+i1;
2387                    const int nz=z+i2;
2388                  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]) {
2389                      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);
2390                      if (i0==0 && i1==0) {                      if (i0==0 && i1==0) {
# Line 2257  void Brick::createPattern() Line 2400  void Brick::createPattern()
2400                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2401                                  if (j>0) {                                  if (j>0) {
2402                                      if (i>0)                                      if (i>0)
2403                                          colIndices[firstDOF+j-1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0, numShared);
2404                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2405                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2406                                          colIndices[firstDOF+j-1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0, numShared);
2407                                  }                                  }
2408                                  if (i>0)                                  if (i>0)
2409                                      colIndices[firstDOF+j-nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0, numShared);
2410                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2411                                  if (i<nDOF1-1)                                  if (i<nDOF1-1)
2412                                      colIndices[firstDOF+j+nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0, numShared);
2413                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2414                                      if (i>0)                                      if (i>0)
2415                                          colIndices[firstDOF+j+1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0, numShared);
2416                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2417                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2418                                          colIndices[firstDOF+j+1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0, numShared);
2419                                  }                                  }
2420                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2421                              }                              }
# Line 2291  void Brick::createPattern() Line 2434  void Brick::createPattern()
2434                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2435                                  if (j>0) {                                  if (j>0) {
2436                                      if (i>0)                                      if (i>0)
2437                                          colIndices[firstDOF+j-1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0*nDOF1, numShared);
2438                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2439                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2440                                          colIndices[firstDOF+j-1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0*nDOF1, numShared);
2441                                  }                                  }
2442                                  if (i>0)                                  if (i>0)
2443                                      colIndices[firstDOF+j-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0*nDOF1, numShared);
2444                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2445                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2446                                      colIndices[firstDOF+j+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0*nDOF1, numShared);
2447                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2448                                      if (i>0)                                      if (i>0)
2449                                          colIndices[firstDOF+j+1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0*nDOF1, numShared);
2450                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2451                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2452                                          colIndices[firstDOF+j+1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0*nDOF1, numShared);
2453                                  }                                  }
2454                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2455                              }                              }
# Line 2325  void Brick::createPattern() Line 2468  void Brick::createPattern()
2468                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2469                                  if (j>0) {                                  if (j>0) {
2470                                      if (i>0)                                      if (i>0)
2471                                          colIndices[firstDOF+(j-1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0-nDOF0*nDOF1, numShared);
2472                                      colIndices[firstDOF+(j-1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0, numShared);
2473                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2474                                          colIndices[firstDOF+(j-1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0+nDOF0*nDOF1, numShared);
2475                                  }                                  }
2476                                  if (i>0)                                  if (i>0)
2477                                      colIndices[firstDOF+j*nDOF0-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0-nDOF0*nDOF1, numShared);
2478                                  colIndices[firstDOF+j*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0, numShared);
2479                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2480                                      colIndices[firstDOF+j*nDOF0+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0+nDOF0*nDOF1, numShared);
2481                                  if (j<nDOF1-1) {                                  if (j<nDOF1-1) {
2482                                      if (i>0)                                      if (i>0)
2483                                          colIndices[firstDOF+(j+1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0-nDOF0*nDOF1, numShared);
2484                                      colIndices[firstDOF+(j+1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0, numShared);
2485                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2486                                          colIndices[firstDOF+(j+1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0+nDOF0*nDOF1, numShared);
2487                                  }                                  }
2488                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2489                              }                              }
# Line 2356  void Brick::createPattern() Line 2499  void Brick::createPattern()
2499                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
2500                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2501                              if (i>0)                              if (i>0)
2502                                  colIndices[firstDOF+i-1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
2503                              colIndices[firstDOF+i].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
2504                              if (i<nDOF0-1)                              if (i<nDOF0-1)
2505                                  colIndices[firstDOF+i+1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
2506                              m_dofMap[firstNode+i]=numDOF+numShared;                              m_dofMap[firstNode+i]=numDOF+numShared;
2507                          }                          }
2508                      } else if (i1==0) {                      } else if (i1==0) {
# Line 2374  void Brick::createPattern() Line 2517  void Brick::createPattern()
2517                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
2518                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2519                              if (i>0)                              if (i>0)
2520                                  colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
2521                              colIndices[firstDOF+i*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
2522                              if (i<nDOF1-1)                              if (i<nDOF1-1)
2523                                  colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
2524                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2525                          }                          }
2526                      } else if (i2==0) {                      } else if (i2==0) {
# Line 2392  void Brick::createPattern() Line 2535  void Brick::createPattern()
2535                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2536                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2537                              if (i>0)                              if (i>0)
2538                                  colIndices[firstDOF+(i-1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0*nDOF1, numShared);
2539                              colIndices[firstDOF+i*nDOF0*nDOF1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0*nDOF1, numShared);
2540                              if (i<nDOF2-1)                              if (i<nDOF2-1)
2541                                  colIndices[firstDOF+(i+1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0*nDOF1, numShared);
2542                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2543                          }                          }
2544                      } else {                      } else {
# Line 2409  void Brick::createPattern() Line 2552  void Brick::createPattern()
2552                          offsetInShared.push_back(offsetInShared.back()+1);                          offsetInShared.push_back(offsetInShared.back()+1);
2553                          sendShared.push_back(dof);                          sendShared.push_back(dof);
2554                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
2555                          colIndices[dof].push_back(numShared);                          doublyLink(colIndices, rowIndices, dof, numShared);
2556                          m_dofMap[node]=numDOF+numShared;                          m_dofMap[node]=numDOF+numShared;
2557                          ++numShared;                          ++numShared;
2558                      }                      }
# Line 2418  void Brick::createPattern() Line 2561  void Brick::createPattern()
2561          }          }
2562      }      }
2563    
2564    #pragma omp parallel for
2565        for (int i = 0; i < numShared; i++) {
2566            std::sort(rowIndices[i].begin(), rowIndices[i].end());
2567        }
2568    
2569      // create connector      // create connector
2570      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
2571              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2572              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2573      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2574              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2575              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2576      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);  
2577    
2578      // useful debug output      // useful debug output
2579      /*      /*
# Line 2502  void Brick::createPattern() Line 2632  void Brick::createPattern()
2632          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2633      }      }
2634      */      */
   
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colPattern);  
     Paso_Pattern_free(rowPattern);  
2635  }  }
2636    
2637  //private  //private
2638  void Brick::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,  void Brick::addToMatrixAndRHS(SystemMatrix* S, escript::Data& F,
2639           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2640           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2641  {  {
# Line 2533  void Brick::addToMatrixAndRHS(Paso_Syste Line 2659  void Brick::addToMatrixAndRHS(Paso_Syste
2659          }          }
2660      }      }
2661      if (addS) {      if (addS) {
2662          addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);          S->add(rowIndex, EM_S);
2663      }      }
2664  }  }
2665    
# Line 2858  void Brick::interpolateNodesOnFaces(escr Line 2984  void Brick::interpolateNodesOnFaces(escr
2984    
2985  namespace  namespace
2986  {  {
2987      // Calculates a guassian blur colvolution matrix for 3D      // Calculates a gaussian blur convolution matrix for 3D
2988      // See wiki article on the subject      // See wiki article on the subject
2989      double* get3DGauss(unsigned radius, double sigma)      double* get3DGauss(unsigned radius, double sigma)
2990      {      {
2991          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)];
2992          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);
2993      double total=0;          double total=0;
2994      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
2995      for (int z=-r;z<=r;++z)          for (int z=-r;z<=r;++z)
2996      {          {
2997          for (int y=-r;y<=r;++y)              for (int y=-r;y<=r;++y)
2998          {              {
2999          for (int x=-r;x<=r;++x)                  for (int x=-r;x<=r;++x)
3000          {                          {        
3001              arr[(x+r)+(y+r)*(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));
3002              total+=arr[(x+r)+(y+r)*(r*2+1)];                      total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];    
3003          }                  }
3004          }              }
3005      }          }
3006      double invtotal=1/total;          double invtotal=1/total;
3007      for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)          for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
3008      {          {
3009          arr[p]*=invtotal;              arr[p]*=invtotal;
3010      }          }
3011      return arr;          return arr;
3012      }      }
3013            
3014      // applies conv to source to get a point.      // applies conv to source to get a point.
# Line 2890  namespace Line 3016  namespace
3016      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)
3017      {      {
3018          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3019      size_t sbase=bx+by*width+bz*width*height;          size_t sbase=bx+by*width+bz*width*height;
3020      double result=0;          double result=0;
3021      for (int z=0;z<2*radius+1;++z)          for (int z=0;z<2*radius+1;++z)
3022      {          {
3023          for (int y=0;y<2*radius+1;++y)              for (int y=0;y<2*radius+1;++y)
3024          {                  {    
3025          for (int x=0;x<2*radius+1;++x)                  for (int x=0;x<2*radius+1;++x)
3026          {                  {
3027              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];
3028          }                  }
3029          }              }
3030      }          }
3031            // use this line for "pass-though" (return the centre point value)
3032    //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3033          return result;                return result;      
3034      }      }
3035  }  }
3036    
3037    /* This is a wrapper for filtered (and non-filtered) randoms
3038     * For detailed doco see randomFillWorker
3039    */
3040    escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
3041           const escript::FunctionSpace& what,
3042           long seed, const boost::python::tuple& filter) const
3043    {
3044        int numvals=escript::DataTypes::noValues(shape);
3045        if (len(filter)>0 && (numvals!=1))
3046        {
3047            throw RipleyException("Ripley only supports filters for scalar data.");
3048        }
3049        escript::Data res=randomFillWorker(shape, seed, filter);
3050        if (res.getFunctionSpace()!=what)
3051        {
3052            escript::Data r=escript::Data(res, what);
3053            return r;
3054        }
3055        return res;
3056    }
3057    
3058  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
3059  The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.  The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.
# Line 2944  inset=2*radius+1 Line 3092  inset=2*radius+1
3092  This is to ensure that values at distance `radius` from the shared/overlapped element  This is to ensure that values at distance `radius` from the shared/overlapped element
3093  that ripley has.  that ripley has.
3094  */  */
3095  escript::Data Brick::randomFill(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
3096  {  {
3097      if (m_numDim!=3)      unsigned int radius=0;  // these are only used by gaussian
3098        double sigma=0.5;
3099        
3100        unsigned int numvals=escript::DataTypes::noValues(shape);
3101        
3102        if (len(filter)==0)
3103      {      {
3104          throw RipleyException("Brick must be 3D.");      // nothing special required here yet
     }  
     if (len(filter)!=3) {  
         throw RipleyException("Unsupported random filter");  
3105      }      }
3106      boost::python::extract<string> ex(filter[0]);      else if (len(filter)==3)
     if (!ex.check() || (ex()!="gaussian"))  
3107      {      {
3108          throw RipleyException("Unsupported random filter");          boost::python::extract<string> ex(filter[0]);
3109            if (!ex.check() || (ex()!="gaussian"))
3110            {
3111                throw RipleyException("Unsupported random filter for Brick.");
3112            }
3113            boost::python::extract<unsigned int> ex1(filter[1]);
3114            if (!ex1.check())
3115            {
3116                throw RipleyException("Radius of gaussian filter must be a positive integer.");
3117            }
3118            radius=ex1();
3119            sigma=0.5;
3120            boost::python::extract<double> ex2(filter[2]);
3121            if (!ex2.check() || (sigma=ex2())<=0)
3122            {
3123                throw RipleyException("Sigma must be a positive floating point number.");
3124            }            
3125      }      }
3126      boost::python::extract<unsigned int> ex1(filter[1]);      else
     if (!ex1.check())  
3127      {      {
3128          throw RipleyException("Radius of gaussian filter must be a positive integer.");          throw RipleyException("Unsupported random filter");
3129      }      }
3130      unsigned int radius=ex1();  
3131      double sigma=0.5;      // number of points in the internal region
3132      boost::python::extract<double> ex2(filter[2]);      // that is, the ones we need smoothed versions of
3133      if (!ex2.check() || (sigma=ex2())<=0)      const dim_t internal[3] = { m_NN[0], m_NN[1], m_NN[2] };
     {  
         throw RipleyException("Sigma must be a postive floating point number.");  
     }      
       
     size_t internal[3];  
     internal[0]=m_NE[0]+1;  // number of points in the internal region  
     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  
3134      size_t ext[3];      size_t ext[3];
3135      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
3136      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3137      ext[2]=internal[2]+2*radius;    // for smoothing      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3138            
3139      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3140      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3141    
3142      if ((2*radius>=internal[0]) || (2*radius>=internal[1]) || (2*radius>=internal[2]))      if (2*radius>=internal[0]-4)
3143      {      {
3144          throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");          throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
3145      }      }
3146            if (2*radius>=internal[1]-4)
3147        {
3148            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3149        }
3150        if (2*radius>=internal[2]-4)
3151        {
3152            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3153        }    
3154        
3155      double* src=new double[ext[0]*ext[1]*ext[2]];      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3156      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3157            
     
3158  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3159          
3160        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3161        {
3162        // since the dimensions are equal for all ranks, this exception
3163        // will be thrown on all ranks
3164        throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3165    
3166        }
3167      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3168      dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];      dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
3169      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3170    #endif    
3171    
3172    /*    
3173        // if we wanted to test a repeating pattern
3174        size_t basex=0;
3175        size_t basey=0;
3176        size_t basez=0;
3177    #ifdef ESYS_MPI    
3178        basex=X*m_gNE[0]/m_NX[0];
3179        basey=Y*m_gNE[1]/m_NX[1];
3180        basez=Z*m_gNE[2]/m_NX[2];
3181        
3182    cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
3183            
3184    #endif    
3185        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3186    */
3187    
3188    #ifdef ESYS_MPI
3189    
3190    
3191    
3192      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);
3193      size_t inset=2*radius+1;          size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3194            // there is an overlap of two points both of which need to have "radius" points on either side.
3195            
3196      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
3197      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;  
3198      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3199            
3200      Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen);          Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);    
3201            
3202      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
3203      MPI_Status stats[50];      MPI_Status stats[50];
3204      short rused=0;      short rused=0;
3205            
# Line 3019  escript::Data Brick::randomFill(long see Line 3210  escript::Data Brick::randomFill(long see
3210      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3211            
3212            
3213      block.copyUsedFromBuffer(src);      block.copyAllToBuffer(src);
       
3214            
3215      int comserr=0;          int comserr=0;    
3216      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
3217      {      {
3218      message& m=incoms[i];          message& m=incoms[i];
3219      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++));
3220      block.setUsed(m.destbuffid);          block.setUsed(m.destbuffid);
3221      }      }
3222    
3223      for (size_t i=0;i<outcoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3224      {      {
3225      message& m=outcoms[i];          message& m=outcoms[i];
3226      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++));
3227      }          }    
3228            
3229      if (!comserr)      if (!comserr)
# Line 3043  escript::Data Brick::randomFill(long see Line 3233  escript::Data Brick::randomFill(long see
3233    
3234      if (comserr)      if (comserr)
3235      {      {
3236      // 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.
3237      // 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.
3238      // 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
3239          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
3240      }      }
3241            
3242      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3243        
3244  #endif      #endif    
3245      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      
3246      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
     // don't need to check for exwrite because we just made it  
     escript::DataVector& dv=resdat.getExpandedVectorReference();  
     double* convolution=get3DGauss(radius, sigma);  
     for (size_t z=0;z<(internal[2]);++z)  
3247      {      {
3248      for (size_t y=0;y<(internal[1]);++y)            
3249      {          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3250          for (size_t x=0;x<(internal[0]);++x)          escript::Data resdat(0, shape, fs , true);
3251          {              // don't need to check for exwrite because we just made it
3252          dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);          escript::DataVector& dv=resdat.getExpandedVectorReference();
3253                
3254          }      
3255      }          // now we need to copy values over
3256            for (size_t z=0;z<(internal[2]);++z)
3257            {
3258                for (size_t y=0;y<(internal[1]);++y)    
3259                {
3260                    for (size_t x=0;x<(internal[0]);++x)
3261                    {
3262                        for (unsigned int i=0;i<numvals;++i)
3263                        {
3264                            dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3265                        }
3266                    }
3267                }
3268            }  
3269            delete[] src;
3270            return resdat;      
3271      }      }
3272      delete[] convolution;      else        // filter enabled  
3273      delete[] src;      {
3274      return resdat;      
3275  }          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3276            escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3277            // don't need to check for exwrite because we just made it
3278            escript::DataVector& dv=resdat.getExpandedVectorReference();
3279            double* convolution=get3DGauss(radius, sigma);
3280    
3281            for (size_t z=0;z<(internal[2]);++z)
3282            {
3283                for (size_t y=0;y<(internal[1]);++y)    
3284                {
3285                    for (size_t x=0;x<(internal[0]);++x)
3286                    {    
3287                        dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3288                
3289                    }
3290                }
3291            }
3292        
3293            delete[] convolution;
3294            delete[] src;
3295            return resdat;
3296        
3297        }
3298    }
3299    
3300  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const
3301    {
3302      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3303      //is the found element even owned by this rank      //is the found element even owned by this rank
3304        // (inside owned or shared elements but will map to an owned element)
3305      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3306          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3307                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
3308            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3309                    + m_dx[dim]/2.;
3310            if (min > coords[dim] || max < coords[dim]) {
3311              return NOT_MINE;              return NOT_MINE;
3312          }          }
3313      }      }
# Line 3091  int Brick::findNode(const double *coords Line 3315  int Brick::findNode(const double *coords
3315      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3316      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3317      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3318        
3319        //check if the point is even inside the domain
3320        if (x < 0 || y < 0 || z < 0
3321                || x > m_length[0] || y > m_length[1] || z > m_length[2])
3322            return NOT_MINE;
3323            
3324      // distance in elements      // distance in elements
3325      int ex = (int) floor(x / m_dx[0]);      int ex = (int) floor(x / m_dx[0]);
3326      int ey = (int) floor(y / m_dx[1]);      int ey = (int) floor(y / m_dx[1]);
# Line 3102  int Brick::findNode(const double *coords Line 3332  int Brick::findNode(const double *coords
3332          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
3333      }      }
3334      //find the closest node      //find the closest node
3335      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
3336          double xdist = x - (ex + dx)*m_dx[0];          double xdist = x - (ex + dx)*m_dx[0];
3337          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
3338              double ydist = y - (ey + dy)*m_dx[1];              double ydist = y - (ey + dy)*m_dx[1];
3339              for (int dz = 0; dz < 2; dz++) {              for (int dz = 0; dz < 1; dz++) {
3340                  double zdist = z - (ez + dz)*m_dx[2];                  double zdist = z - (ez + dz)*m_dx[2];
3341                  double total = xdist*xdist + ydist*ydist + zdist*zdist;                  double total = xdist*xdist + ydist*ydist + zdist*zdist;
3342                  if (total < minDist) {                  if (total < minDist) {
3343                      closest = INDEX3(ex+dy, ey+dy, ez+dz, m_NE[0]+1, m_NE[1]+1);                      closest = INDEX3(ex+dy-m_offset[0], ey+dy-m_offset[1],
3344                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3345                      minDist = total;                      minDist = total;
3346                  }                  }
3347              }              }
3348          }          }
3349      }      }
3350        if (closest == NOT_MINE) {
3351            throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");
3352        }
3353      return closest;      return closest;
3354  }  }
3355    
3356  void Brick::setAssembler(std::string type, std::map<std::string,  Assembler_ptr Brick::createAssembler(std::string type, std::map<std::string,
3357          escript::Data> constants) {          escript::Data> constants) const
3358      if (type.compare("WaveAssembler") == 0) {  {
3359          delete assembler;      if (type.compare("DefaultAssembler") == 0) {
3360          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));
3361        } else if (type.compare("WaveAssembler") == 0) {
3362            return Assembler_ptr(new WaveAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN, constants));
3363        } else if (type.compare("LameAssembler") == 0) {
3364            return Assembler_ptr(new LameAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
3365      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3366          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Brick does not support the"
3367                                  " requested assembler");                                  " requested assembler");
3368      }      }
3369  }  }

Legend:
Removed from v.4653  
changed lines
  Added in v.4988

  ViewVC Help
Powered by ViewVC 1.1.26