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

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

  ViewVC Help
Powered by ViewVC 1.1.26