/[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

revision 4650 by jfenwick, Wed Feb 5 04:16:01 2014 UTC revision 4861 by sshaw, Thu Apr 10 05:17:47 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>
20  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
21  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
22  #include <ripley/DefaultAssembler3D.h>  #include <ripley/DefaultAssembler3D.h>
23  #include <ripley/WaveAssembler3D.h>  #include <ripley/WaveAssembler3D.h>
24    #include <ripley/LameAssembler3D.h>
25    #include <ripley/domainhelpers.h>
26  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
27    
28  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 42  using esysUtils::FileWriter; Line 47  using esysUtils::FileWriter;
47    
48  namespace ripley {  namespace ripley {
49    
50    int indexOfMax(int a, int b, int c) {
51        if (a > b) {
52            if (c > a) {
53                return 2;
54            }
55            return 0;
56        } else if (b > c) {
57            return 1;
58        }
59        return 2;
60    }
61    
62  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,
63               double x1, double y1, double z1, int d0, int d1, int d2,               double x1, double y1, double z1, int d0, int d1, int d2,
64               const std::vector<double>& points, const std::vector<int>& tags,               const std::vector<double>& points, const std::vector<int>& tags,
65               const simap_t& tagnamestonums) :               const simap_t& tagnamestonums) :
66      RipleyDomain(3)      RipleyDomain(3)
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 240  Brick::Brick(int n0, int n1, int n2, dou Line 214  Brick::Brick(int n0, int n1, int n2, dou
214    
215  Brick::~Brick()  Brick::~Brick()
216  {  {
     Paso_SystemMatrixPattern_free(m_pattern);  
     Paso_Connector_free(m_connector);  
217      delete assembler;      delete assembler;
218  }  }
219    
# Line 411  void Brick::readNcGrid(escript::Data& ou Line 383  void Brick::readNcGrid(escript::Data& ou
383  #endif  #endif
384  }  }
385    
386    #ifdef USE_BOOSTIO
387    void Brick::readBinaryGridFromZipped(escript::Data& out, string filename,
388                               const ReaderParameters& params) const
389    {
390        // the mapping is not universally correct but should work on our
391        // supported platforms
392        switch (params.dataType) {
393            case DATATYPE_INT32:
394                readBinaryGridZippedImpl<int>(out, filename, params);
395                break;
396            case DATATYPE_FLOAT32:
397                readBinaryGridZippedImpl<float>(out, filename, params);
398                break;
399            case DATATYPE_FLOAT64:
400                readBinaryGridZippedImpl<double>(out, filename, params);
401                break;
402            default:
403                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
404        }
405    }
406    #endif
407    
408  void Brick::readBinaryGrid(escript::Data& out, string filename,  void Brick::readBinaryGrid(escript::Data& out, string filename,
409                             const ReaderParameters& params) const                             const ReaderParameters& params) const
410  {  {
# Line 547  void Brick::readBinaryGridImpl(escript:: Line 541  void Brick::readBinaryGridImpl(escript::
541      f.close();      f.close();
542  }  }
543    
544    #ifdef USE_BOOSTIO
545    template<typename ValueType>
546    void Brick::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
547                                   const ReaderParameters& params) const
548    {
549        // check destination function space
550        int myN0, myN1, myN2;
551        if (out.getFunctionSpace().getTypeCode() == Nodes) {
552            myN0 = m_NN[0];
553            myN1 = m_NN[1];
554            myN2 = m_NN[2];
555        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
556                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
557            myN0 = m_NE[0];
558            myN1 = m_NE[1];
559            myN2 = m_NE[2];
560        } else
561            throw RipleyException("readBinaryGridFromZipped(): invalid function space for output data object");
562    
563        if (params.first.size() != 3)
564            throw RipleyException("readBinaryGridFromZipped(): argument 'first' must have 3 entries");
565    
566        if (params.numValues.size() != 3)
567            throw RipleyException("readBinaryGridFromZipped(): argument 'numValues' must have 3 entries");
568    
569        if (params.multiplier.size() != 3)
570            throw RipleyException("readBinaryGridFromZipped(): argument 'multiplier' must have 3 entries");
571        for (size_t i=0; i<params.multiplier.size(); i++)
572            if (params.multiplier[i]<1)
573                throw RipleyException("readBinaryGridFromZipped(): all multipliers must be positive");
574    
575        // check file existence and size
576        ifstream f(filename.c_str(), ifstream::binary);
577        if (f.fail()) {
578            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
579        }
580        f.seekg(0, ios::end);
581        const int numComp = out.getDataPointSize();
582        int filesize = f.tellg();
583        f.seekg(0, ios::beg);
584        std::vector<char> compressed(filesize);
585        f.read((char*)&compressed[0], filesize);
586        f.close();
587        std::vector<char> decompressed = unzip(compressed);
588        filesize = decompressed.size();
589        const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
590        if (filesize < reqsize) {
591            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
592        }
593    
594        // check if this rank contributes anything
595        if (params.first[0] >= m_offset[0]+myN0 ||
596                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
597                params.first[1] >= m_offset[1]+myN1 ||
598                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
599                params.first[2] >= m_offset[2]+myN2 ||
600                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
601            return;
602        }
603    
604        // now determine how much this rank has to write
605    
606        // first coordinates in data object to write to
607        const int first0 = max(0, params.first[0]-m_offset[0]);
608        const int first1 = max(0, params.first[1]-m_offset[1]);
609        const int first2 = max(0, params.first[2]-m_offset[2]);
610        // indices to first value in file
611        const int idx0 = max(0, m_offset[0]-params.first[0]);
612        const int idx1 = max(0, m_offset[1]-params.first[1]);
613        const int idx2 = max(0, m_offset[2]-params.first[2]);
614        // number of values to read
615        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
616        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
617        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
618    
619        out.requireWrite();
620        vector<ValueType> values(num0*numComp);
621        const int dpp = out.getNumDataPointsPerSample();
622    
623        for (int z=0; z<num2; z++) {
624            for (int y=0; y<num1; y++) {
625                const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
626                                 +(idx2+z)*params.numValues[0]*params.numValues[1]);
627                memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
628                
629                for (int x=0; x<num0; x++) {
630                    const int baseIndex = first0+x*params.multiplier[0]
631                                         +(first1+y*params.multiplier[1])*myN0
632                                         +(first2+z*params.multiplier[2])*myN0*myN1;
633                    for (int m2=0; m2<params.multiplier[2]; m2++) {
634                        for (int m1=0; m1<params.multiplier[1]; m1++) {
635                            for (int m0=0; m0<params.multiplier[0]; m0++) {
636                                const int dataIndex = baseIndex+m0
637                                               +m1*myN0
638                                               +m2*myN0*myN1;
639                                double* dest = out.getSampleDataRW(dataIndex);
640                                for (int c=0; c<numComp; c++) {
641                                    ValueType val = values[x*numComp+c];
642    
643                                    if (params.byteOrder != BYTEORDER_NATIVE) {
644                                        char* cval = reinterpret_cast<char*>(&val);
645                                        // this will alter val!!
646                                        byte_swap32(cval);
647                                    }
648                                    if (!std::isnan(val)) {
649                                        for (int q=0; q<dpp; q++) {
650                                            *dest++ = static_cast<double>(val);
651                                        }
652                                    }
653                                }
654                            }
655                        }
656                    }
657                }
658            }
659        }
660    }
661    #endif
662    
663  void Brick::writeBinaryGrid(const escript::Data& in, string filename,  void Brick::writeBinaryGrid(const escript::Data& in, string filename,
664                              int byteOrder, int dataType) const                              int byteOrder, int dataType) const
665  {  {
# Line 786  const int* Brick::borrowSampleReferenceI Line 899  const int* Brick::borrowSampleReferenceI
899          case FaceElements:          case FaceElements:
900          case ReducedFaceElements:          case ReducedFaceElements:
901              return &m_faceId[0];              return &m_faceId[0];
902            case Points:
903                return &m_diracPointNodeIDs[0];
904          default:          default:
905              break;              break;
906      }      }
# Line 1944  void Brick::nodesToDOF(escript::Data& ou Line 2059  void Brick::nodesToDOF(escript::Data& ou
2059  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const
2060  {  {
2061      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
2062      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2063      // 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
2064      const_cast<escript::Data*>(&in)->expand();      const_cast<escript::Data*>(&in)->expand();
2065      Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));      coupler->startCollect(in.getSampleDataRO(0));
2066    
2067      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2068      out.requireWrite();      out.requireWrite();
2069      const double* buffer = Paso_Coupler_finishCollect(coupler);      const double* buffer = coupler->finishCollect();
2070    
2071  #pragma omp parallel for  #pragma omp parallel for
2072      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 2075  void Brick::dofToNodes(escript::Data& ou
2075                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
2076          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
2077      }      }
     Paso_Coupler_free(coupler);  
2078  }  }
2079    
2080  //private  //private
# Line 1981  void Brick::populateSampleIds() Line 2095  void Brick::populateSampleIds()
2095          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2096      }      }
2097      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2098      m_nodeId.resize(getNumNodes());      
2099      m_dofId.resize(numDOF);      try {
2100      m_elementId.resize(getNumElements());          m_nodeId.resize(getNumNodes());
2101            m_dofId.resize(numDOF);
2102            m_elementId.resize(getNumElements());
2103        } catch (const std::length_error& le) {
2104            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2105        }
2106        
2107      // populate face element counts      // populate face element counts
2108      //left      //left
2109      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2228  void Brick::createPattern() Line 2347  void Brick::createPattern()
2347      RankVector neighbour;      RankVector neighbour;
2348      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
2349      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
2350      int numShared=0;      int numShared=0, expectedShared=0;;
2351      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
2352      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];
2353      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 2361  void Brick::createPattern()
2361                  const int nx=x+i0;                  const int nx=x+i0;
2362                  const int ny=y+i1;                  const int ny=y+i1;
2363                  const int nz=z+i2;                  const int nz=z+i2;
2364                    if (!(nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2])) {
2365                        continue;
2366                    }
2367                    if (i0==0 && i1==0)
2368                        expectedShared += nDOF0*nDOF1;
2369                    else if (i0==0 && i2==0)
2370                        expectedShared += nDOF0*nDOF2;
2371                    else if (i1==0 && i2==0)
2372                        expectedShared += nDOF1*nDOF2;
2373                    else if (i0==0)
2374                        expectedShared += nDOF0;
2375                    else if (i1==0)
2376                        expectedShared += nDOF1;
2377                    else if (i2==0)
2378                        expectedShared += nDOF2;
2379                    else
2380                        expectedShared++;
2381                }
2382            }
2383        }
2384        
2385        vector<IndexVector> rowIndices(expectedShared);
2386        
2387        for (int i2=-1; i2<2; i2++) {
2388            for (int i1=-1; i1<2; i1++) {
2389                for (int i0=-1; i0<2; i0++) {
2390                    // skip this rank
2391                    if (i0==0 && i1==0 && i2==0)
2392                        continue;
2393                    // location of neighbour rank
2394                    const int nx=x+i0;
2395                    const int ny=y+i1;
2396                    const int nz=z+i2;
2397                  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]) {
2398                      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);
2399                      if (i0==0 && i1==0) {                      if (i0==0 && i1==0) {
# Line 2257  void Brick::createPattern() Line 2409  void Brick::createPattern()
2409                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2410                                  if (j>0) {                                  if (j>0) {
2411                                      if (i>0)                                      if (i>0)
2412                                          colIndices[firstDOF+j-1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0, numShared);
2413                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2414                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2415                                          colIndices[firstDOF+j-1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0, numShared);
2416                                  }                                  }
2417                                  if (i>0)                                  if (i>0)
2418                                      colIndices[firstDOF+j-nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0, numShared);
2419                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2420                                  if (i<nDOF1-1)                                  if (i<nDOF1-1)
2421                                      colIndices[firstDOF+j+nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0, numShared);
2422                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2423                                      if (i>0)                                      if (i>0)
2424                                          colIndices[firstDOF+j+1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0, numShared);
2425                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2426                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2427                                          colIndices[firstDOF+j+1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0, numShared);
2428                                  }                                  }
2429                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2430                              }                              }
# Line 2291  void Brick::createPattern() Line 2443  void Brick::createPattern()
2443                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2444                                  if (j>0) {                                  if (j>0) {
2445                                      if (i>0)                                      if (i>0)
2446                                          colIndices[firstDOF+j-1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0*nDOF1, numShared);
2447                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2448                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2449                                          colIndices[firstDOF+j-1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0*nDOF1, numShared);
2450                                  }                                  }
2451                                  if (i>0)                                  if (i>0)
2452                                      colIndices[firstDOF+j-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0*nDOF1, numShared);
2453                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2454                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2455                                      colIndices[firstDOF+j+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0*nDOF1, numShared);
2456                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2457                                      if (i>0)                                      if (i>0)
2458                                          colIndices[firstDOF+j+1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0*nDOF1, numShared);
2459                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2460                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2461                                          colIndices[firstDOF+j+1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0*nDOF1, numShared);
2462                                  }                                  }
2463                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2464                              }                              }
# Line 2325  void Brick::createPattern() Line 2477  void Brick::createPattern()
2477                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2478                                  if (j>0) {                                  if (j>0) {
2479                                      if (i>0)                                      if (i>0)
2480                                          colIndices[firstDOF+(j-1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0-nDOF0*nDOF1, numShared);
2481                                      colIndices[firstDOF+(j-1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0, numShared);
2482                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2483                                          colIndices[firstDOF+(j-1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0+nDOF0*nDOF1, numShared);
2484                                  }                                  }
2485                                  if (i>0)                                  if (i>0)
2486                                      colIndices[firstDOF+j*nDOF0-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0-nDOF0*nDOF1, numShared);
2487                                  colIndices[firstDOF+j*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0, numShared);
2488                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2489                                      colIndices[firstDOF+j*nDOF0+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0+nDOF0*nDOF1, numShared);
2490                                  if (j<nDOF1-1) {                                  if (j<nDOF1-1) {
2491                                      if (i>0)                                      if (i>0)
2492                                          colIndices[firstDOF+(j+1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0-nDOF0*nDOF1, numShared);
2493                                      colIndices[firstDOF+(j+1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0, numShared);
2494                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2495                                          colIndices[firstDOF+(j+1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0+nDOF0*nDOF1, numShared);
2496                                  }                                  }
2497                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2498                              }                              }
# Line 2356  void Brick::createPattern() Line 2508  void Brick::createPattern()
2508                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
2509                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2510                              if (i>0)                              if (i>0)
2511                                  colIndices[firstDOF+i-1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
2512                              colIndices[firstDOF+i].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
2513                              if (i<nDOF0-1)                              if (i<nDOF0-1)
2514                                  colIndices[firstDOF+i+1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
2515                              m_dofMap[firstNode+i]=numDOF+numShared;                              m_dofMap[firstNode+i]=numDOF+numShared;
2516                          }                          }
2517                      } else if (i1==0) {                      } else if (i1==0) {
# Line 2374  void Brick::createPattern() Line 2526  void Brick::createPattern()
2526                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
2527                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2528                              if (i>0)                              if (i>0)
2529                                  colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
2530                              colIndices[firstDOF+i*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
2531                              if (i<nDOF1-1)                              if (i<nDOF1-1)
2532                                  colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
2533                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2534                          }                          }
2535                      } else if (i2==0) {                      } else if (i2==0) {
# Line 2392  void Brick::createPattern() Line 2544  void Brick::createPattern()
2544                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2545                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2546                              if (i>0)                              if (i>0)
2547                                  colIndices[firstDOF+(i-1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0*nDOF1, numShared);
2548                              colIndices[firstDOF+i*nDOF0*nDOF1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0*nDOF1, numShared);
2549                              if (i<nDOF2-1)                              if (i<nDOF2-1)
2550                                  colIndices[firstDOF+(i+1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0*nDOF1, numShared);
2551                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2552                          }                          }
2553                      } else {                      } else {
# Line 2409  void Brick::createPattern() Line 2561  void Brick::createPattern()
2561                          offsetInShared.push_back(offsetInShared.back()+1);                          offsetInShared.push_back(offsetInShared.back()+1);
2562                          sendShared.push_back(dof);                          sendShared.push_back(dof);
2563                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
2564                          colIndices[dof].push_back(numShared);                          doublyLink(colIndices, rowIndices, dof, numShared);
2565                          m_dofMap[node]=numDOF+numShared;                          m_dofMap[node]=numDOF+numShared;
2566                          ++numShared;                          ++numShared;
2567                      }                      }
# Line 2418  void Brick::createPattern() Line 2570  void Brick::createPattern()
2570          }          }
2571      }      }
2572    
2573    #pragma omp parallel for
2574        for (int i = 0; i < numShared; i++) {
2575            std::sort(rowIndices[i].begin(), rowIndices[i].end());
2576        }
2577    
2578      // create connector      // create connector
2579      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
2580              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2581              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2582      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2583              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2584              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2585      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);  
2586    
2587      // create main and couple blocks      // create main and couple blocks
2588      Paso_Pattern *mainPattern = createMainPattern();      paso::Pattern_ptr mainPattern = createMainPattern();
2589      Paso_Pattern *colPattern, *rowPattern;      paso::Pattern_ptr colPattern, rowPattern;
2590      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, colPattern, rowPattern);
2591    
2592      // allocate paso distribution      // allocate paso distribution
2593      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2594              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2595    
2596      // finally create the system matrix      // finally create the system matrix
2597      m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,      m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2598              distribution, distribution, mainPattern, colPattern, rowPattern,              distribution, distribution, mainPattern, colPattern, rowPattern,
2599              m_connector, m_connector);              m_connector, m_connector));
   
     Paso_Distribution_free(distribution);  
2600    
2601      // useful debug output      // useful debug output
2602      /*      /*
# Line 2502  void Brick::createPattern() Line 2655  void Brick::createPattern()
2655          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2656      }      }
2657      */      */
   
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colPattern);  
     Paso_Pattern_free(rowPattern);  
2658  }  }
2659    
2660  //private  //private
2661  void Brick::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,  void Brick::addToMatrixAndRHS(paso::SystemMatrix_ptr S, escript::Data& F,
2662           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2663           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2664  {  {
# Line 2858  void Brick::interpolateNodesOnFaces(escr Line 3007  void Brick::interpolateNodesOnFaces(escr
3007    
3008  namespace  namespace
3009  {  {
3010      // Calculates a guassian blur colvolution matrix for 3D      // Calculates a gaussian blur convolution matrix for 3D
3011      // See wiki article on the subject      // See wiki article on the subject
3012      double* get3DGauss(unsigned radius, double sigma)      double* get3DGauss(unsigned radius, double sigma)
3013      {      {
3014          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)];
3015          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);
3016      double total=0;          double total=0;
3017      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
3018      for (int z=-r;z<=r;++z)          for (int z=-r;z<=r;++z)
3019      {          {
3020          for (int y=-r;y<=r;++y)              for (int y=-r;y<=r;++y)
3021          {              {
3022          for (int x=-r;x<=r;++x)                  for (int x=-r;x<=r;++x)
3023          {                          {        
3024              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));
3025              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)];    
3026          }                  }
3027          }              }
3028      }          }
3029      double invtotal=1/total;          double invtotal=1/total;
3030      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)
3031      {          {
3032          arr[p]*=invtotal;              arr[p]*=invtotal;
3033      }          }
3034      return arr;          return arr;
3035      }      }
3036            
3037      // applies conv to source to get a point.      // applies conv to source to get a point.
# Line 2890  namespace Line 3039  namespace
3039      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)
3040      {      {
3041          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3042      size_t sbase=bx+by*width+bz*width*height;          size_t sbase=bx+by*width+bz*width*height;
3043      double result=0;          double result=0;
3044      for (int z=0;z<2*radius+1;++z)          for (int z=0;z<2*radius+1;++z)
3045      {          {
3046          for (int y=0;y<2*radius+1;++y)              for (int y=0;y<2*radius+1;++y)
3047          {                  {    
3048          for (int x=0;x<2*radius+1;++x)                  for (int x=0;x<2*radius+1;++x)
3049          {                  {
3050              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];
3051          }                  }
3052          }              }
3053      }          }
3054            // use this line for "pass-though" (return the centre point value)
3055    //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3056          return result;                return result;      
3057      }      }
3058  }  }
3059    
3060    /* This is a wrapper for filtered (and non-filtered) randoms
3061     * For detailed doco see randomFillWorker
3062    */
3063    escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
3064           const escript::FunctionSpace& what,
3065           long seed, const boost::python::tuple& filter) const
3066    {
3067        int numvals=escript::DataTypes::noValues(shape);
3068        if (len(filter)>0 && (numvals!=1))
3069        {
3070            throw RipleyException("Ripley only supports filters for scalar data.");
3071        }
3072        escript::Data res=randomFillWorker(shape, seed, filter);
3073        if (res.getFunctionSpace()!=what)
3074        {
3075            escript::Data r=escript::Data(res, what);
3076            return r;
3077        }
3078        return res;
3079    }
3080    
3081  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
3082  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 3115  inset=2*radius+1
3115  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
3116  that ripley has.  that ripley has.
3117  */  */
3118  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
3119  {  {
3120      if (m_numDim!=3)      unsigned int radius=0;  // these are only used by gaussian
3121        double sigma=0.5;
3122        
3123        unsigned int numvals=escript::DataTypes::noValues(shape);
3124        
3125        if (len(filter)==0)
3126      {      {
3127          throw RipleyException("Brick must be 3D.");      // nothing special required here yet
     }  
     if (len(filter)!=3) {  
         throw RipleyException("Unsupported random filter");  
3128      }      }
3129      boost::python::extract<string> ex(filter[0]);      else if (len(filter)==3)
     if (!ex.check() || (ex()!="gaussian"))  
3130      {      {
3131          throw RipleyException("Unsupported random filter");          boost::python::extract<string> ex(filter[0]);
3132            if (!ex.check() || (ex()!="gaussian"))
3133            {
3134                throw RipleyException("Unsupported random filter for Brick.");
3135            }
3136            boost::python::extract<unsigned int> ex1(filter[1]);
3137            if (!ex1.check())
3138            {
3139                throw RipleyException("Radius of gaussian filter must be a positive integer.");
3140            }
3141            radius=ex1();
3142            sigma=0.5;
3143            boost::python::extract<double> ex2(filter[2]);
3144            if (!ex2.check() || (sigma=ex2())<=0)
3145            {
3146                throw RipleyException("Sigma must be a positive floating point number.");
3147            }            
3148      }      }
3149      boost::python::extract<unsigned int> ex1(filter[1]);      else
     if (!ex1.check())  
3150      {      {
3151          throw RipleyException("Radius of gaussian filter must be a positive integer.");          throw RipleyException("Unsupported random filter");
3152      }      }
3153      unsigned int radius=ex1();  
3154      double sigma=0.5;      // number of points in the internal region
3155      boost::python::extract<double> ex2(filter[2]);      // that is, the ones we need smoothed versions of
3156      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  
3157      size_t ext[3];      size_t ext[3];
3158      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
3159      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3160      ext[2]=internal[2]+2*radius;    // for smoothing      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3161            
3162      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3163      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3164    
3165      if ((2*radius>=internal[0]) || (2*radius>=internal[1]) || (2*radius>=internal[2]))      if (2*radius>=internal[0]-4)
3166      {      {
3167          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");
3168      }      }
3169            if (2*radius>=internal[1]-4)
3170        {
3171            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3172        }
3173        if (2*radius>=internal[2]-4)
3174        {
3175            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3176        }    
3177        
3178      double* src=new double[ext[0]*ext[1]*ext[2]];      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3179      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3180            
     
3181  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3182          
3183        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3184        {
3185        // since the dimensions are equal for all ranks, this exception
3186        // will be thrown on all ranks
3187        throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3188    
3189        }
3190      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3191      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];
3192      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3193    #endif    
3194    
3195    /*    
3196        // if we wanted to test a repeating pattern
3197        size_t basex=0;
3198        size_t basey=0;
3199        size_t basez=0;
3200    #ifdef ESYS_MPI    
3201        basex=X*m_gNE[0]/m_NX[0];
3202        basey=Y*m_gNE[1]/m_NX[1];
3203        basez=Z*m_gNE[2]/m_NX[2];
3204            
3205    cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
3206        
3207    #endif    
3208        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3209    */
3210    
3211    #ifdef ESYS_MPI
3212    
3213    
3214    
3215      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);
3216      size_t inset=2*radius+1;          size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3217            // there is an overlap of two points both of which need to have "radius" points on either side.
3218            
3219      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
3220      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;  
3221      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3222            
3223      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);    
3224            
3225      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
3226      MPI_Status stats[50];      MPI_Status stats[50];
3227      short rused=0;      short rused=0;
3228            
# Line 3019  escript::Data Brick::randomFill(long see Line 3233  escript::Data Brick::randomFill(long see
3233      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3234            
3235            
3236      block.copyUsedFromBuffer(src);      block.copyAllToBuffer(src);
       
3237            
3238      int comserr=0;          int comserr=0;    
3239      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
3240      {      {
3241      message& m=incoms[i];          message& m=incoms[i];
3242      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++));
3243      block.setUsed(m.destbuffid);          block.setUsed(m.destbuffid);
3244      }      }
3245    
3246      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3247      {      {
3248      message& m=incoms[i];          message& m=outcoms[i];
3249      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++));
3250      }          }    
3251            
3252      if (!comserr)      if (!comserr)
# Line 3043  escript::Data Brick::randomFill(long see Line 3256  escript::Data Brick::randomFill(long see
3256    
3257      if (comserr)      if (comserr)
3258      {      {
3259      // 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.
3260      // 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.
3261      // 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
3262          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
3263      }      }
3264            
3265      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3266        
3267  #endif      #endif    
3268      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      
3269      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
3270      // don't need to check for exwrite because we just made it      {
3271      escript::DataVector& dv=resdat.getExpandedVectorReference();        
3272      double* convolution=get3DGauss(radius, sigma);          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3273      for (size_t z=0;z<(internal[2]);++z)          escript::Data resdat(0, shape, fs , true);
3274            // don't need to check for exwrite because we just made it
3275            escript::DataVector& dv=resdat.getExpandedVectorReference();
3276        
3277        
3278            // now we need to copy values over
3279            for (size_t z=0;z<(internal[2]);++z)
3280            {
3281                for (size_t y=0;y<(internal[1]);++y)    
3282                {
3283                    for (size_t x=0;x<(internal[0]);++x)
3284                    {
3285                        for (unsigned int i=0;i<numvals;++i)
3286                        {
3287                            dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3288                        }
3289                    }
3290                }
3291            }  
3292            delete[] src;
3293            return resdat;      
3294        }
3295        else        // filter enabled  
3296      {      {
3297      for (size_t y=0;y<(internal[1]);++y)          
3298      {          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3299          for (size_t x=0;x<(internal[0]);++x)          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3300          {              // don't need to check for exwrite because we just made it
3301          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();
3302                    double* convolution=get3DGauss(radius, sigma);
3303          }  
3304      }          for (size_t z=0;z<(internal[2]);++z)
3305      }          {
3306      delete[] convolution;              for (size_t y=0;y<(internal[1]);++y)    
3307      delete[] src;              {
3308      return resdat;                  for (size_t x=0;x<(internal[0]);++x)
3309                    {    
3310                        dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3311                
3312                    }
3313                }
3314            }
3315        
3316            delete[] convolution;
3317            delete[] src;
3318            return resdat;
3319        
3320        }
3321  }  }
3322    
3323    
# Line 3081  escript::Data Brick::randomFill(long see Line 3328  escript::Data Brick::randomFill(long see
3328  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const {
3329      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3330      //is the found element even owned by this rank      //is the found element even owned by this rank
3331        // (inside owned or shared elements but will map to an owned element)
3332      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3333          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3334                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
3335            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3336                    + m_dx[dim]/2.;
3337            if (min > coords[dim] || max < coords[dim]) {
3338              return NOT_MINE;              return NOT_MINE;
3339          }          }
3340      }      }
# Line 3091  int Brick::findNode(const double *coords Line 3342  int Brick::findNode(const double *coords
3342      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3343      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3344      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3345        
3346        //check if the point is even inside the domain
3347        if (x < 0 || y < 0 || z < 0
3348                || x > m_length[0] || y > m_length[1] || z > m_length[2])
3349            return NOT_MINE;
3350            
3351      // distance in elements      // distance in elements
3352      int ex = (int) floor(x / m_dx[0]);      int ex = (int) floor(x / m_dx[0]);
3353      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 3359  int Brick::findNode(const double *coords
3359          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
3360      }      }
3361      //find the closest node      //find the closest node
3362      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
3363          double xdist = x - (ex + dx)*m_dx[0];          double xdist = x - (ex + dx)*m_dx[0];
3364          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
3365              double ydist = y - (ey + dy)*m_dx[1];              double ydist = y - (ey + dy)*m_dx[1];
3366              for (int dz = 0; dz < 2; dz++) {              for (int dz = 0; dz < 1; dz++) {
3367                  double zdist = z - (ez + dz)*m_dx[2];                  double zdist = z - (ez + dz)*m_dx[2];
3368                  double total = xdist*xdist + ydist*ydist + zdist*zdist;                  double total = xdist*xdist + ydist*ydist + zdist*zdist;
3369                  if (total < minDist) {                  if (total < minDist) {
3370                      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],
3371                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3372                      minDist = total;                      minDist = total;
3373                  }                  }
3374              }              }
3375          }          }
3376      }      }
3377        if (closest == NOT_MINE) {
3378            throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");
3379        }
3380      return closest;      return closest;
3381  }  }
3382    
3383  void Brick::setAssembler(std::string type, std::map<std::string,  void Brick::setAssembler(std::string type, std::map<std::string,
3384          escript::Data> constants) {          escript::Data> constants) {
3385      if (type.compare("WaveAssembler") == 0) {      if (type.compare("WaveAssembler") == 0) {
3386            if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3387                throw RipleyException("Domain already using a different custom assembler");
3388            assembler_type = WAVE_ASSEMBLER;
3389          delete assembler;          delete assembler;
3390          assembler = new WaveAssembler3D(this, m_dx, m_NX, m_NE, m_NN, constants);          assembler = new WaveAssembler3D(this, m_dx, m_NX, m_NE, m_NN, constants);
3391        } else if (type.compare("LameAssembler") == 0) {
3392            if (assembler_type != LAME_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3393                throw RipleyException("Domain already using a different custom assembler");
3394            assembler_type = LAME_ASSEMBLER;
3395            delete assembler;
3396            assembler = new LameAssembler3D(this, m_dx, m_NX, m_NE, m_NN);
3397      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3398          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Brick does not support the"
3399                                  " requested assembler");                                  " requested assembler");
3400      }      }
3401  }  }

Legend:
Removed from v.4650  
changed lines
  Added in v.4861

  ViewVC Help
Powered by ViewVC 1.1.26