/[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 4629 by sshaw, Fri Jan 24 03:29:25 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 33  Line 37 
37    
38  #include <iomanip>  #include <iomanip>
39    
40    #include "esysUtils/EsysRandom.h"
41    #include "blocktools.h"
42    
43    
44  using namespace std;  using namespace std;
45  using esysUtils::FileWriter;  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 std::map<std::string, int>& 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 51  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 225  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 236  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 407  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 456  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 488  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 505  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 539  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 782  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 1879  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 1940  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      Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));      // expand data object if necessary to be able to grab the whole data
2055        const_cast<escript::Data*>(&in)->expand();
2056        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 1954  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 1975  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 2222  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 2236  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 2251  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 2285  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 2319  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 2350  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 2368  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 2386  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 2403  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 2412  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 2496  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 2527  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 2850  void Brick::interpolateNodesOnFaces(escr Line 2982  void Brick::interpolateNodesOnFaces(escr
2982      }      }
2983  }  }
2984    
2985  int Brick::findNode(const double *coords) const {  namespace
2986    {
2987        // Calculates a gaussian blur convolution matrix for 3D
2988        // See wiki article on the subject
2989        double* get3DGauss(unsigned radius, double sigma)
2990        {
2991            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);
2993            double total=0;
2994            int r=static_cast<int>(radius);
2995            for (int z=-r;z<=r;++z)
2996            {
2997                for (int y=-r;y<=r;++y)
2998                {
2999                    for (int x=-r;x<=r;++x)
3000                    {        
3001                        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)+(z+r)*(r*2+1)*(r*2+1)];    
3003                    }
3004                }
3005            }
3006            double invtotal=1/total;
3007            for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
3008            {
3009                arr[p]*=invtotal;
3010            }
3011            return arr;
3012        }
3013        
3014        // applies conv to source to get a point.
3015        // (xp, yp) are the coords in the source matrix not the destination matrix
3016        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;
3019            size_t sbase=bx+by*width+bz*width*height;
3020            double result=0;
3021            for (int z=0;z<2*radius+1;++z)
3022            {
3023                for (int y=0;y<2*radius+1;++y)
3024                {    
3025                    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];
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;      
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.
3059    The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.
3060    A parameter radius  gives the size of the stencil used for the smoothing.
3061    A point on the left hand edge for example, will still require `radius` extra points to the left
3062    in order to complete the stencil.
3063    
3064    All local calculation is done on an array called `src`, with
3065    dimensions = ext[0] * ext[1] *ext[2].
3066    Where ext[i]= internal[i]+2*radius.
3067    
3068    Now for MPI there is overlap to deal with. We need to share both the overlapping
3069    values themselves but also the external region.
3070    
3071    In a hypothetical 1-D case:
3072    
3073    
3074    1234567
3075    would be split into two ranks thus:
3076    123(4)  (4)567     [4 being a shared element]
3077    
3078    If the radius is 2.   There will be padding elements on the outside:
3079    
3080    pp123(4)  (4)567pp
3081    
3082    To ensure that 4 can be correctly computed on both ranks, values from the other rank
3083    need to be known.
3084    
3085    pp123(4)56   23(4)567pp
3086    
3087    Now in our case, we wout set all the values 23456 on the left rank and send them to the
3088    right hand rank.
3089    
3090    So the edges _may_ need to be shared at a distance `inset` from all boundaries.
3091    inset=2*radius+1    
3092    This is to ensure that values at distance `radius` from the shared/overlapped element
3093    that ripley has.
3094    */
3095    escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
3096    {
3097        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        // nothing special required here yet
3105        }
3106        else if (len(filter)==3)
3107        {
3108            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        else
3127        {
3128            throw RipleyException("Unsupported random filter");
3129        }
3130    
3131        // number of points in the internal region
3132        // that is, the ones we need smoothed versions of
3133        const dim_t internal[3] = { m_NN[0], m_NN[1], m_NN[2] };
3134        size_t ext[3];
3135        ext[0]=(size_t)internal[0]+2*radius;  // includes points we need as input
3136        ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3137        ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3138        
3139        // now we check to see if the radius is acceptable
3140        // That is, would not cross multiple ranks in MPI
3141    
3142        if (2*radius>=internal[0]-4)
3143        {
3144            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]*numvals];
3156        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3157        
3158    #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];
3168        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]);
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);
3193        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
3197        size_t ymidlen=ext[1]-2*inset;  
3198        size_t zmidlen=ext[2]-2*inset;
3199        
3200        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
3203        MPI_Status stats[50];
3204        short rused=0;
3205        
3206        messvec incoms;
3207        messvec outcoms;
3208        
3209        grid.generateInNeighbours(X, Y, Z ,incoms);
3210        grid.generateOutNeighbours(X, Y, Z, outcoms);
3211        
3212        
3213        block.copyAllToBuffer(src);
3214        
3215        int comserr=0;    
3216        for (size_t i=0;i<incoms.size();++i)
3217        {
3218            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++));
3220            block.setUsed(m.destbuffid);
3221        }
3222    
3223        for (size_t i=0;i<outcoms.size();++i)
3224        {
3225            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++));
3227        }    
3228        
3229        if (!comserr)
3230        {    
3231            comserr=MPI_Waitall(rused, reqs, stats);
3232        }
3233    
3234        if (comserr)
3235        {
3236        // 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.
3238        // however, we have no reason to believe coms work at this point anyway
3239            throw RipleyException("Error in coms for randomFill");      
3240        }
3241        
3242        block.copyUsedFromBuffer(src);
3243        
3244    #endif    
3245        
3246        if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
3247        {
3248          
3249            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3250            escript::Data resdat(0, shape, fs , true);
3251            // don't need to check for exwrite because we just made it
3252            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        else        // filter enabled  
3273        {
3274        
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
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 2863  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 2874  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.4629  
changed lines
  Added in v.4988

  ViewVC Help
Powered by ViewVC 1.1.26