/[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 4687 by jfenwick, Wed Feb 19 00:03:29 2014 UTC revision 4861 by sshaw, Thu Apr 10 05:17:47 2014 UTC
# Line 14  Line 14 
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 43  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 56  Brick::Brick(int n0, int n1, int n2, dou Line 81  Brick::Brick(int n0, int n1, int n2, dou
81          d2=1;          d2=1;
82      }      }
83      bool warn=false;      bool warn=false;
84      // if number of subdivisions is non-positive, try to subdivide by the same  
85      // ratio as the number of elements      std::vector<int> factors;
86      if (d0<=0 && d1<=0 && d2<=0) {      int ranks = m_mpiInfo->size;
87          warn=true;      int epr[3] = {n0,n1,n2};
88          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));      int d[3] = {d0,d1,d2};
89          d0=max(1, d0);      if (d0<=0 || d1<=0 || d2<=0) {
90          d1=max(1, (int)(d0*n1/(float)n0));          for (int i = 0; i < 3; i++) {
91          d2=m_mpiInfo->size/(d0*d1);              if (d[i] < 1) {
92          if (d0*d1*d2 != m_mpiInfo->size) {                  d[i] = 1;
93              // ratios not the same so leave "smallest" side undivided and try                  continue;
94              // dividing 2 sides only              }
95              if (n0>=n1) {              epr[i] = -1; // can no longer be max
96                  if (n1>=n2) {              if (ranks % d[i] != 0) {
97                      d0=d1=0;                  throw RipleyException("Invalid number of spatial subdivisions");
98                      d2=1;              }
99                  } else {              //remove
100                      d0=d2=0;              ranks /= d[i];
101                      d1=1;          }
102                  }          factorise(factors, ranks);
103              } else {          if (factors.size() != 0) {
104                  if (n0>=n2) {              warn = true;
105                      d0=d1=0;          }
106                      d2=1;      }
107                  } else {      while (factors.size() > 0) {
108                      d0=1;          int i = indexOfMax(epr[0],epr[1],epr[2]);
109                      d1=d2=0;          int f = factors.back();
110                  }          factors.pop_back();
111              }          d[i] *= f;
112          }          epr[i] /= f;
     }  
     if (d0<=0 && d1<=0) {  
         warn=true;  
         d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));  
         d1=m_mpiInfo->size/d0;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n0>n1) {  
                 d0=0;  
                 d1=1;  
             } else {  
                 d0=1;  
                 d1=0;  
             }  
         }  
     } else if (d0<=0 && d2<=0) {  
         warn=true;  
         d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1))));  
         d2=m_mpiInfo->size/d0;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n0>n2) {  
                 d0=0;  
                 d2=1;  
             } else {  
                 d0=1;  
                 d2=0;  
             }  
         }  
     } else if (d1<=0 && d2<=0) {  
         warn=true;  
         d1=max(1, int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1))));  
         d2=m_mpiInfo->size/d1;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n1>n2) {  
                 d1=0;  
                 d2=1;  
             } else {  
                 d1=1;  
                 d2=0;  
             }  
         }  
     }  
     if (d0<=0) {  
         // d1,d2 are preset, determine d0  
         d0=m_mpiInfo->size/(d1*d2);  
     } else if (d1<=0) {  
         // d0,d2 are preset, determine d1  
         d1=m_mpiInfo->size/(d0*d2);  
     } else if (d2<=0) {  
         // d0,d1 are preset, determine d2  
         d2=m_mpiInfo->size/(d0*d1);  
113      }      }
114        d0 = d[0]; d1 = d[1]; d2 = d[2];
115    
116      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
117      // among number of ranks      // among number of ranks
118      if (d0*d1*d2 != m_mpiInfo->size)      if (d0*d1*d2 != m_mpiInfo->size){
119          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
120        }
121      if (warn) {      if (warn) {
122          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
123              << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;              << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;
# Line 241  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 412  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 548  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 1947  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 1963  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 1984  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 2231  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 2245  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 2260  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 2294  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 2328  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 2359  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 2377  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 2395  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 2412  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 2421  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 2505  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 2861  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)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));                      arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));
3025              total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];                      total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];    
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)*(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 2893  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)          // use this line for "pass-though" (return the centre point value)
3055  //  return source[sbase+(radius)+(radius)*width+(radius)*width*height];  //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3056      return result;                return result;      
3057      }      }
3058  }  }
3059    
# Line 2921  escript::Data Brick::randomFill(const es Line 3067  escript::Data Brick::randomFill(const es
3067      int numvals=escript::DataTypes::noValues(shape);      int numvals=escript::DataTypes::noValues(shape);
3068      if (len(filter)>0 && (numvals!=1))      if (len(filter)>0 && (numvals!=1))
3069      {      {
3070      throw RipleyException("Ripley only supports filters for scalar data.");          throw RipleyException("Ripley only supports filters for scalar data.");
3071      }      }
3072      escript::Data res=randomFillWorker(shape, seed, filter);      escript::Data res=randomFillWorker(shape, seed, filter);
3073      if (res.getFunctionSpace()!=what)      if (res.getFunctionSpace()!=what)
3074      {      {
3075      escript::Data r=escript::Data(res, what);          escript::Data r=escript::Data(res, what);
3076      return r;          return r;
3077      }      }
3078      return res;      return res;
3079  }  }
# Line 2971  that ripley has. Line 3117  that ripley has.
3117  */  */
3118  escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const  escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
3119  {  {
3120      if (m_numDim!=3)      unsigned int radius=0;  // these are only used by gaussian
     {  
         throw RipleyException("Brick must be 3D.");  
     }  
       
     unsigned int radius=0;  // these are only used by gaussian  
3121      double sigma=0.5;      double sigma=0.5;
3122            
3123      unsigned int numvals=escript::DataTypes::noValues(shape);      unsigned int numvals=escript::DataTypes::noValues(shape);
3124            
3125      if (len(filter)==0)      if (len(filter)==0)
3126      {      {
3127      // nothing special required here yet      // nothing special required here yet
3128      }      }
3129      else if (len(filter)==3)      else if (len(filter)==3)
3130      {      {
3131      boost::python::extract<string> ex(filter[0]);          boost::python::extract<string> ex(filter[0]);
3132      if (!ex.check() || (ex()!="gaussian"))          if (!ex.check() || (ex()!="gaussian"))
3133      {          {
3134          throw RipleyException("Unsupported random filter for Brick.");              throw RipleyException("Unsupported random filter for Brick.");
3135      }          }
3136      boost::python::extract<unsigned int> ex1(filter[1]);          boost::python::extract<unsigned int> ex1(filter[1]);
3137      if (!ex1.check())          if (!ex1.check())
3138      {          {
3139          throw RipleyException("Radius of gaussian filter must be a positive integer.");              throw RipleyException("Radius of gaussian filter must be a positive integer.");
3140      }          }
3141      radius=ex1();          radius=ex1();
3142      sigma=0.5;          sigma=0.5;
3143      boost::python::extract<double> ex2(filter[2]);          boost::python::extract<double> ex2(filter[2]);
3144      if (!ex2.check() || (sigma=ex2())<=0)          if (!ex2.check() || (sigma=ex2())<=0)
3145      {          {
3146          throw RipleyException("Sigma must be a postive floating point number.");              throw RipleyException("Sigma must be a positive floating point number.");
3147      }                      }            
3148      }      }
3149      else      else
3150      {      {
3151          throw RipleyException("Unsupported random filter");          throw RipleyException("Unsupported random filter");
3152      }      }
       
       
3153    
3154            // number of points in the internal region
3155      size_t internal[3];      // that is, the ones we need smoothed versions of
3156      internal[0]=m_NE[0]+1;  // number of points in the internal region      const dim_t internal[3] = { m_NN[0], m_NN[1], m_NN[2] };
     internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of  
     internal[2]=m_NE[2]+1;  // that is, the ones we need smoothed versions of  
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
# Line 3045  escript::Data Brick::randomFillWorker(co Line 3182  escript::Data Brick::randomFillWorker(co
3182        
3183      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3184      {      {
3185      // since the dimensions are equal for all ranks, this exception      // since the dimensions are equal for all ranks, this exception
3186      // will be thrown on all ranks      // will be thrown on all ranks
3187      throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");      throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3188    
3189      }      }
3190      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
# Line 3063  escript::Data Brick::randomFillWorker(co Line 3200  escript::Data Brick::randomFillWorker(co
3200  #ifdef ESYS_MPI      #ifdef ESYS_MPI    
3201      basex=X*m_gNE[0]/m_NX[0];      basex=X*m_gNE[0]/m_NX[0];
3202      basey=Y*m_gNE[1]/m_NX[1];      basey=Y*m_gNE[1]/m_NX[1];
3203      basez=Z*m_gNE[2]/m_NX[2];          basez=Z*m_gNE[2]/m_NX[2];
 #endif      
     if (seed==0)  
     {  
     seed=2; // since we are using the seed parameter as the spacing and 0 spacing causes an exception  
     }  
     esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, seed, basex, basey, basez);  
 */  
       
3204            
3205  /*  cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
 cout << "Pattern:\n";      
 for (int i=0;i<ext[0]*ext[1]*ext[2];)  
 {  
     cout << src[i++] << " ";  
     if (i%ext[0]==0)  
       cout << "\n";  
     if (i%(ext[0]*ext[1])==0)  
       cout << "\n";  
 }*/  
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  #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+2;    // Its +2 not +1 because a whole element is shared (and hence      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3217          // there is an overlap of two points both of which need to have "radius" points on either side.          // there is an overlap of two points both of which need to have "radius" points on either side.
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, numvals);          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 3112  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3235  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3235            
3236      block.copyAllToBuffer(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<outcoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3247      {      {
3248      message& m=outcoms[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 3134  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3256  for (int i=0;i<ext[0]*ext[1]*ext[2];)
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            
3269  /*          if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
 cout << "Pattern (post transfer):\n";      
 for (int i=0;i<ext[0]*ext[1]*ext[2];)  
 {  
     cout << src[i++] << " ";  
     if (i%ext[0]==0)  
       cout << "\n";  
     if (i%(ext[0]*ext[1])==0)  
       cout << "\n";  
 }   */  
       
     if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe  
3270      {      {
3271                
3272      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3273      escript::Data resdat(0, shape, fs , true);          escript::Data resdat(0, shape, fs , true);
3274      // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3275      escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3276            
3277            
3278      // now we need to copy values over          // now we need to copy values over
3279      for (size_t z=0;z<(internal[2]);++z)          for (size_t z=0;z<(internal[2]);++z)
3280      {          {
3281          for (size_t y=0;y<(internal[1]);++y)                  for (size_t y=0;y<(internal[1]);++y)    
3282          {              {
3283          for (size_t x=0;x<(internal[0]);++x)                  for (size_t x=0;x<(internal[0]);++x)
3284          {                  {
3285              for (unsigned int i=0;i<numvals;++i)                      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];                          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;          delete[] src;
3293      return resdat;                return resdat;      
3294      }      }
3295      else        // filter enabled        else        // filter enabled  
3296      {      {
3297            
3298      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3299      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3300      // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3301      escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3302      double* convolution=get3DGauss(radius, sigma);          double* convolution=get3DGauss(radius, sigma);
3303    
3304  // cout << "Convolution matrix\n";          for (size_t z=0;z<(internal[2]);++z)
3305  // size_t di=(radius*2+1);          {
3306  // double ctot=0;              for (size_t y=0;y<(internal[1]);++y)    
3307  // for (int i=0;i<di*di*di;++i)              {
3308  // {                  for (size_t x=0;x<(internal[0]);++x)
3309  //     cout << convolution[i] << " ";                  {    
3310  //     ctot+=convolution[i];                      dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3311  //     if ((i+1)%di==0)              
3312  //     {                  }
3313  //  cout << "\n";              }
3314  //     }          }
3315  //     if ((i+1)%(di*di)==0)      
3316  //     {          delete[] convolution;
3317  //  cout << "\n";          delete[] src;
3318  //     }          return resdat;
 // }  
 //  
 // cout << "Sum of matrix is = " << ctot << endl;  
   
     for (size_t z=0;z<(internal[2]);++z)  
     {  
         for (size_t y=0;y<(internal[1]);++y)      
         {  
         for (size_t x=0;x<(internal[0]);++x)  
         {      
             dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);  
               
         }  
         }  
     }  
       
 // cout << "\nResulting matrix:\n";  
 //     for (size_t z=0;z<(internal[2]);++z)  
 //     {  
 //  for (size_t y=0;y<(internal[1]);++y)      
 //  {  
 //      for (size_t x=0;x<(internal[0]);++x)  
 //      {      
 //      cout << dv[x+y*(internal[0])+z*internal[0]*internal[1]] << " ";  
 //      }  
 //      cout << endl;  
 //  }  
 //  cout << endl;  
 //     }  
   
       
     delete[] convolution;  
     delete[] src;  
     return resdat;  
3319            
3320      }      }
3321  }  }
# Line 3268  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 3303  int Brick::findNode(const double *coords Line 3383  int Brick::findNode(const double *coords
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.4687  
changed lines
  Added in v.4861

  ViewVC Help
Powered by ViewVC 1.1.26