/[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 4817 by caltinay, Fri Mar 28 08:04:09 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    
# Line 18  Line 19 
19  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
20  #include <ripley/DefaultAssembler3D.h>  #include <ripley/DefaultAssembler3D.h>
21  #include <ripley/WaveAssembler3D.h>  #include <ripley/WaveAssembler3D.h>
22    #include <ripley/LameAssembler3D.h>
23    #include <ripley/domainhelpers.h>
24  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
25    
26  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 42  using esysUtils::FileWriter; Line 45  using esysUtils::FileWriter;
45    
46  namespace ripley {  namespace ripley {
47    
48    int indexOfMax(int a, int b, int c) {
49        if (a > b) {
50            if (c > a) {
51                return 2;
52            }
53            return 0;
54        } else if (b > c) {
55            return 1;
56        }
57        return 2;
58    }
59    
60  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,
61               double x1, double y1, double z1, int d0, int d1, int d2,               double x1, double y1, double z1, int d0, int d1, int d2,
62               const std::vector<double>& points, const std::vector<int>& tags,               const std::vector<double>& points, const std::vector<int>& tags,
# Line 55  Brick::Brick(int n0, int n1, int n2, dou Line 70  Brick::Brick(int n0, int n1, int n2, dou
70          d2=1;          d2=1;
71      }      }
72      bool warn=false;      bool warn=false;
73      // if number of subdivisions is non-positive, try to subdivide by the same  
74      // ratio as the number of elements      std::vector<int> factors;
75      if (d0<=0 && d1<=0 && d2<=0) {      int ranks = m_mpiInfo->size;
76          warn=true;      int epr[3] = {n0,n1,n2};
77          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));      int d[3] = {d0,d1,d2};
78          d0=max(1, d0);      if (d0<=0 || d1<=0 || d2<=0) {
79          d1=max(1, (int)(d0*n1/(float)n0));          for (int i = 0; i < 3; i++) {
80          d2=m_mpiInfo->size/(d0*d1);              if (d[i] < 1) {
81          if (d0*d1*d2 != m_mpiInfo->size) {                  d[i] = 1;
82              // ratios not the same so leave "smallest" side undivided and try                  continue;
83              // dividing 2 sides only              }
84              if (n0>=n1) {              epr[i] = -1; // can no longer be max
85                  if (n1>=n2) {              if (ranks % d[i] != 0) {
86                      d0=d1=0;                  throw RipleyException("Invalid number of spatial subdivisions");
87                      d2=1;              }
88                  } else {              //remove
89                      d0=d2=0;              ranks /= d[i];
90                      d1=1;          }
91                  }          factorise(factors, ranks);
92              } else {          if (factors.size() != 0) {
93                  if (n0>=n2) {              warn = true;
94                      d0=d1=0;          }
95                      d2=1;      }
96                  } else {      while (factors.size() > 0) {
97                      d0=1;          int i = indexOfMax(epr[0],epr[1],epr[2]);
98                      d1=d2=0;          int f = factors.back();
99                  }          factors.pop_back();
100              }          d[i] *= f;
101          }          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);  
102      }      }
103        d0 = d[0]; d1 = d[1]; d2 = d[2];
104    
105      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
106      // among number of ranks      // among number of ranks
107      if (d0*d1*d2 != m_mpiInfo->size)      if (d0*d1*d2 != m_mpiInfo->size){
108          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
109        }
110      if (warn) {      if (warn) {
111          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
112              << 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 203  Brick::Brick(int n0, int n1, int n2, dou
203    
204  Brick::~Brick()  Brick::~Brick()
205  {  {
206      Paso_SystemMatrixPattern_free(m_pattern);      paso::SystemMatrixPattern_free(m_pattern);
     Paso_Connector_free(m_connector);  
207      delete assembler;      delete assembler;
208  }  }
209    
# Line 411  void Brick::readNcGrid(escript::Data& ou Line 373  void Brick::readNcGrid(escript::Data& ou
373  #endif  #endif
374  }  }
375    
376    #ifdef USE_BOOSTIO
377    void Brick::readBinaryGridFromZipped(escript::Data& out, string filename,
378                               const ReaderParameters& params) const
379    {
380        // the mapping is not universally correct but should work on our
381        // supported platforms
382        switch (params.dataType) {
383            case DATATYPE_INT32:
384                readBinaryGridZippedImpl<int>(out, filename, params);
385                break;
386            case DATATYPE_FLOAT32:
387                readBinaryGridZippedImpl<float>(out, filename, params);
388                break;
389            case DATATYPE_FLOAT64:
390                readBinaryGridZippedImpl<double>(out, filename, params);
391                break;
392            default:
393                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
394        }
395    }
396    #endif
397    
398  void Brick::readBinaryGrid(escript::Data& out, string filename,  void Brick::readBinaryGrid(escript::Data& out, string filename,
399                             const ReaderParameters& params) const                             const ReaderParameters& params) const
400  {  {
# Line 547  void Brick::readBinaryGridImpl(escript:: Line 531  void Brick::readBinaryGridImpl(escript::
531      f.close();      f.close();
532  }  }
533    
534    #ifdef USE_BOOSTIO
535    template<typename ValueType>
536    void Brick::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
537                                   const ReaderParameters& params) const
538    {
539        // check destination function space
540        int myN0, myN1, myN2;
541        if (out.getFunctionSpace().getTypeCode() == Nodes) {
542            myN0 = m_NN[0];
543            myN1 = m_NN[1];
544            myN2 = m_NN[2];
545        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
546                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
547            myN0 = m_NE[0];
548            myN1 = m_NE[1];
549            myN2 = m_NE[2];
550        } else
551            throw RipleyException("readBinaryGridFromZipped(): invalid function space for output data object");
552    
553        if (params.first.size() != 3)
554            throw RipleyException("readBinaryGridFromZipped(): argument 'first' must have 3 entries");
555    
556        if (params.numValues.size() != 3)
557            throw RipleyException("readBinaryGridFromZipped(): argument 'numValues' must have 3 entries");
558    
559        if (params.multiplier.size() != 3)
560            throw RipleyException("readBinaryGridFromZipped(): argument 'multiplier' must have 3 entries");
561        for (size_t i=0; i<params.multiplier.size(); i++)
562            if (params.multiplier[i]<1)
563                throw RipleyException("readBinaryGridFromZipped(): all multipliers must be positive");
564    
565        // check file existence and size
566        ifstream f(filename.c_str(), ifstream::binary);
567        if (f.fail()) {
568            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
569        }
570        f.seekg(0, ios::end);
571        const int numComp = out.getDataPointSize();
572        int filesize = f.tellg();
573        f.seekg(0, ios::beg);
574        std::vector<char> compressed(filesize);
575        f.read((char*)&compressed[0], filesize);
576        f.close();
577        std::vector<char> decompressed = unzip(compressed);
578        filesize = decompressed.size();
579        const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
580        if (filesize < reqsize) {
581            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
582        }
583    
584        // check if this rank contributes anything
585        if (params.first[0] >= m_offset[0]+myN0 ||
586                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
587                params.first[1] >= m_offset[1]+myN1 ||
588                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
589                params.first[2] >= m_offset[2]+myN2 ||
590                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
591            return;
592        }
593    
594        // now determine how much this rank has to write
595    
596        // first coordinates in data object to write to
597        const int first0 = max(0, params.first[0]-m_offset[0]);
598        const int first1 = max(0, params.first[1]-m_offset[1]);
599        const int first2 = max(0, params.first[2]-m_offset[2]);
600        // indices to first value in file
601        const int idx0 = max(0, m_offset[0]-params.first[0]);
602        const int idx1 = max(0, m_offset[1]-params.first[1]);
603        const int idx2 = max(0, m_offset[2]-params.first[2]);
604        // number of values to read
605        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
606        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
607        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
608    
609        out.requireWrite();
610        vector<ValueType> values(num0*numComp);
611        const int dpp = out.getNumDataPointsPerSample();
612    
613        for (int z=0; z<num2; z++) {
614            for (int y=0; y<num1; y++) {
615                const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
616                                 +(idx2+z)*params.numValues[0]*params.numValues[1]);
617                memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
618                
619                for (int x=0; x<num0; x++) {
620                    const int baseIndex = first0+x*params.multiplier[0]
621                                         +(first1+y*params.multiplier[1])*myN0
622                                         +(first2+z*params.multiplier[2])*myN0*myN1;
623                    for (int m2=0; m2<params.multiplier[2]; m2++) {
624                        for (int m1=0; m1<params.multiplier[1]; m1++) {
625                            for (int m0=0; m0<params.multiplier[0]; m0++) {
626                                const int dataIndex = baseIndex+m0
627                                               +m1*myN0
628                                               +m2*myN0*myN1;
629                                double* dest = out.getSampleDataRW(dataIndex);
630                                for (int c=0; c<numComp; c++) {
631                                    ValueType val = values[x*numComp+c];
632    
633                                    if (params.byteOrder != BYTEORDER_NATIVE) {
634                                        char* cval = reinterpret_cast<char*>(&val);
635                                        // this will alter val!!
636                                        byte_swap32(cval);
637                                    }
638                                    if (!std::isnan(val)) {
639                                        for (int q=0; q<dpp; q++) {
640                                            *dest++ = static_cast<double>(val);
641                                        }
642                                    }
643                                }
644                            }
645                        }
646                    }
647                }
648            }
649        }
650    }
651    #endif
652    
653  void Brick::writeBinaryGrid(const escript::Data& in, string filename,  void Brick::writeBinaryGrid(const escript::Data& in, string filename,
654                              int byteOrder, int dataType) const                              int byteOrder, int dataType) const
655  {  {
# Line 786  const int* Brick::borrowSampleReferenceI Line 889  const int* Brick::borrowSampleReferenceI
889          case FaceElements:          case FaceElements:
890          case ReducedFaceElements:          case ReducedFaceElements:
891              return &m_faceId[0];              return &m_faceId[0];
892            case Points:
893                return &m_diracPointNodeIDs[0];
894          default:          default:
895              break;              break;
896      }      }
# Line 1944  void Brick::nodesToDOF(escript::Data& ou Line 2049  void Brick::nodesToDOF(escript::Data& ou
2049  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const
2050  {  {
2051      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
2052      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2053      // 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
2054      const_cast<escript::Data*>(&in)->expand();      const_cast<escript::Data*>(&in)->expand();
2055      Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));      coupler->startCollect(in.getSampleDataRO(0));
2056    
2057      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2058      out.requireWrite();      out.requireWrite();
2059      const double* buffer = Paso_Coupler_finishCollect(coupler);      const double* buffer = coupler->finishCollect();
2060    
2061  #pragma omp parallel for  #pragma omp parallel for
2062      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 2065  void Brick::dofToNodes(escript::Data& ou
2065                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
2066          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
2067      }      }
     Paso_Coupler_free(coupler);  
2068  }  }
2069    
2070  //private  //private
# Line 2228  void Brick::createPattern() Line 2332  void Brick::createPattern()
2332      RankVector neighbour;      RankVector neighbour;
2333      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
2334      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
2335      int numShared=0;      int numShared=0, expectedShared=0;;
2336      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
2337      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];
2338      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 2346  void Brick::createPattern()
2346                  const int nx=x+i0;                  const int nx=x+i0;
2347                  const int ny=y+i1;                  const int ny=y+i1;
2348                  const int nz=z+i2;                  const int nz=z+i2;
2349                    if (!(nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2])) {
2350                        continue;
2351                    }
2352                    if (i0==0 && i1==0)
2353                        expectedShared += nDOF0*nDOF1;
2354                    else if (i0==0 && i2==0)
2355                        expectedShared += nDOF0*nDOF2;
2356                    else if (i1==0 && i2==0)
2357                        expectedShared += nDOF1*nDOF2;
2358                    else if (i0==0)
2359                        expectedShared += nDOF0;
2360                    else if (i1==0)
2361                        expectedShared += nDOF1;
2362                    else if (i2==0)
2363                        expectedShared += nDOF2;
2364                    else
2365                        expectedShared++;
2366                }
2367            }
2368        }
2369        
2370        vector<IndexVector> rowIndices(expectedShared);
2371        
2372        for (int i2=-1; i2<2; i2++) {
2373            for (int i1=-1; i1<2; i1++) {
2374                for (int i0=-1; i0<2; i0++) {
2375                    // skip this rank
2376                    if (i0==0 && i1==0 && i2==0)
2377                        continue;
2378                    // location of neighbour rank
2379                    const int nx=x+i0;
2380                    const int ny=y+i1;
2381                    const int nz=z+i2;
2382                  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]) {
2383                      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);
2384                      if (i0==0 && i1==0) {                      if (i0==0 && i1==0) {
# Line 2257  void Brick::createPattern() Line 2394  void Brick::createPattern()
2394                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2395                                  if (j>0) {                                  if (j>0) {
2396                                      if (i>0)                                      if (i>0)
2397                                          colIndices[firstDOF+j-1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0, numShared);
2398                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2399                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2400                                          colIndices[firstDOF+j-1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0, numShared);
2401                                  }                                  }
2402                                  if (i>0)                                  if (i>0)
2403                                      colIndices[firstDOF+j-nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0, numShared);
2404                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2405                                  if (i<nDOF1-1)                                  if (i<nDOF1-1)
2406                                      colIndices[firstDOF+j+nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0, numShared);
2407                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2408                                      if (i>0)                                      if (i>0)
2409                                          colIndices[firstDOF+j+1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0, numShared);
2410                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2411                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2412                                          colIndices[firstDOF+j+1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0, numShared);
2413                                  }                                  }
2414                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2415                              }                              }
# Line 2291  void Brick::createPattern() Line 2428  void Brick::createPattern()
2428                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2429                                  if (j>0) {                                  if (j>0) {
2430                                      if (i>0)                                      if (i>0)
2431                                          colIndices[firstDOF+j-1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0*nDOF1, numShared);
2432                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2433                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2434                                          colIndices[firstDOF+j-1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0*nDOF1, numShared);
2435                                  }                                  }
2436                                  if (i>0)                                  if (i>0)
2437                                      colIndices[firstDOF+j-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0*nDOF1, numShared);
2438                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2439                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2440                                      colIndices[firstDOF+j+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0*nDOF1, numShared);
2441                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2442                                      if (i>0)                                      if (i>0)
2443                                          colIndices[firstDOF+j+1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0*nDOF1, numShared);
2444                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2445                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2446                                          colIndices[firstDOF+j+1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0*nDOF1, numShared);
2447                                  }                                  }
2448                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2449                              }                              }
# Line 2325  void Brick::createPattern() Line 2462  void Brick::createPattern()
2462                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2463                                  if (j>0) {                                  if (j>0) {
2464                                      if (i>0)                                      if (i>0)
2465                                          colIndices[firstDOF+(j-1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0-nDOF0*nDOF1, numShared);
2466                                      colIndices[firstDOF+(j-1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0, numShared);
2467                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2468                                          colIndices[firstDOF+(j-1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0+nDOF0*nDOF1, numShared);
2469                                  }                                  }
2470                                  if (i>0)                                  if (i>0)
2471                                      colIndices[firstDOF+j*nDOF0-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0-nDOF0*nDOF1, numShared);
2472                                  colIndices[firstDOF+j*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0, numShared);
2473                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2474                                      colIndices[firstDOF+j*nDOF0+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0+nDOF0*nDOF1, numShared);
2475                                  if (j<nDOF1-1) {                                  if (j<nDOF1-1) {
2476                                      if (i>0)                                      if (i>0)
2477                                          colIndices[firstDOF+(j+1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0-nDOF0*nDOF1, numShared);
2478                                      colIndices[firstDOF+(j+1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0, numShared);
2479                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2480                                          colIndices[firstDOF+(j+1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0+nDOF0*nDOF1, numShared);
2481                                  }                                  }
2482                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2483                              }                              }
# Line 2356  void Brick::createPattern() Line 2493  void Brick::createPattern()
2493                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
2494                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2495                              if (i>0)                              if (i>0)
2496                                  colIndices[firstDOF+i-1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
2497                              colIndices[firstDOF+i].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
2498                              if (i<nDOF0-1)                              if (i<nDOF0-1)
2499                                  colIndices[firstDOF+i+1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
2500                              m_dofMap[firstNode+i]=numDOF+numShared;                              m_dofMap[firstNode+i]=numDOF+numShared;
2501                          }                          }
2502                      } else if (i1==0) {                      } else if (i1==0) {
# Line 2374  void Brick::createPattern() Line 2511  void Brick::createPattern()
2511                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
2512                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2513                              if (i>0)                              if (i>0)
2514                                  colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
2515                              colIndices[firstDOF+i*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
2516                              if (i<nDOF1-1)                              if (i<nDOF1-1)
2517                                  colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
2518                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2519                          }                          }
2520                      } else if (i2==0) {                      } else if (i2==0) {
# Line 2392  void Brick::createPattern() Line 2529  void Brick::createPattern()
2529                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2530                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2531                              if (i>0)                              if (i>0)
2532                                  colIndices[firstDOF+(i-1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0*nDOF1, numShared);
2533                              colIndices[firstDOF+i*nDOF0*nDOF1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0*nDOF1, numShared);
2534                              if (i<nDOF2-1)                              if (i<nDOF2-1)
2535                                  colIndices[firstDOF+(i+1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0*nDOF1, numShared);
2536                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2537                          }                          }
2538                      } else {                      } else {
# Line 2409  void Brick::createPattern() Line 2546  void Brick::createPattern()
2546                          offsetInShared.push_back(offsetInShared.back()+1);                          offsetInShared.push_back(offsetInShared.back()+1);
2547                          sendShared.push_back(dof);                          sendShared.push_back(dof);
2548                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
2549                          colIndices[dof].push_back(numShared);                          doublyLink(colIndices, rowIndices, dof, numShared);
2550                          m_dofMap[node]=numDOF+numShared;                          m_dofMap[node]=numDOF+numShared;
2551                          ++numShared;                          ++numShared;
2552                      }                      }
# Line 2418  void Brick::createPattern() Line 2555  void Brick::createPattern()
2555          }          }
2556      }      }
2557    
2558    #pragma omp parallel for
2559        for (int i = 0; i < numShared; i++) {
2560            std::sort(rowIndices[i].begin(), rowIndices[i].end());
2561        }
2562    
2563      // create connector      // create connector
2564      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
2565              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2566              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2567      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2568              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2569              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2570      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);  
2571    
2572      // create main and couple blocks      // create main and couple blocks
2573      Paso_Pattern *mainPattern = createMainPattern();      paso::Pattern *mainPattern = createMainPattern();
2574      Paso_Pattern *colPattern, *rowPattern;      paso::Pattern *colPattern, *rowPattern;
2575      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, &colPattern, &rowPattern);
2576    
2577      // allocate paso distribution      // allocate paso distribution
2578      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2579              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2580    
2581      // finally create the system matrix      // finally create the system matrix
2582      m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,      m_pattern = new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2583              distribution, distribution, mainPattern, colPattern, rowPattern,              distribution, distribution, mainPattern, colPattern, rowPattern,
2584              m_connector, m_connector);              m_connector, m_connector);
2585    
     Paso_Distribution_free(distribution);  
   
2586      // useful debug output      // useful debug output
2587      /*      /*
2588      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
# Line 2503  void Brick::createPattern() Line 2641  void Brick::createPattern()
2641      }      }
2642      */      */
2643    
2644      Paso_Pattern_free(mainPattern);      paso::Pattern_free(mainPattern);
2645      Paso_Pattern_free(colPattern);      paso::Pattern_free(colPattern);
2646      Paso_Pattern_free(rowPattern);      paso::Pattern_free(rowPattern);
2647  }  }
2648    
2649  //private  //private
# Line 2858  void Brick::interpolateNodesOnFaces(escr Line 2996  void Brick::interpolateNodesOnFaces(escr
2996    
2997  namespace  namespace
2998  {  {
2999      // Calculates a guassian blur colvolution matrix for 3D      // Calculates a gaussian blur convolution matrix for 3D
3000      // See wiki article on the subject      // See wiki article on the subject
3001      double* get3DGauss(unsigned radius, double sigma)      double* get3DGauss(unsigned radius, double sigma)
3002      {      {
3003          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)];
3004          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);
3005      double total=0;          double total=0;
3006      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
3007      for (int z=-r;z<=r;++z)          for (int z=-r;z<=r;++z)
3008      {          {
3009          for (int y=-r;y<=r;++y)              for (int y=-r;y<=r;++y)
3010          {              {
3011          for (int x=-r;x<=r;++x)                  for (int x=-r;x<=r;++x)
3012          {                          {        
3013              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));
3014              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)];    
3015          }                  }
3016          }              }
3017      }          }
3018      double invtotal=1/total;          double invtotal=1/total;
3019      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)
3020      {          {
3021          arr[p]*=invtotal;              arr[p]*=invtotal;
3022      }          }
3023      return arr;          return arr;
3024      }      }
3025            
3026      // applies conv to source to get a point.      // applies conv to source to get a point.
# Line 2890  namespace Line 3028  namespace
3028      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)
3029      {      {
3030          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3031      size_t sbase=bx+by*width+bz*width*height;          size_t sbase=bx+by*width+bz*width*height;
3032      double result=0;          double result=0;
3033      for (int z=0;z<2*radius+1;++z)          for (int z=0;z<2*radius+1;++z)
3034      {          {
3035          for (int y=0;y<2*radius+1;++y)              for (int y=0;y<2*radius+1;++y)
3036          {                  {    
3037          for (int x=0;x<2*radius+1;++x)                  for (int x=0;x<2*radius+1;++x)
3038          {                  {
3039              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];
3040          }                  }
3041          }              }
3042      }          }
3043            // use this line for "pass-though" (return the centre point value)
3044    //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3045          return result;                return result;      
3046      }      }
3047  }  }
3048    
3049    /* This is a wrapper for filtered (and non-filtered) randoms
3050     * For detailed doco see randomFillWorker
3051    */
3052    escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
3053           const escript::FunctionSpace& what,
3054           long seed, const boost::python::tuple& filter) const
3055    {
3056        int numvals=escript::DataTypes::noValues(shape);
3057        if (len(filter)>0 && (numvals!=1))
3058        {
3059            throw RipleyException("Ripley only supports filters for scalar data.");
3060        }
3061        escript::Data res=randomFillWorker(shape, seed, filter);
3062        if (res.getFunctionSpace()!=what)
3063        {
3064            escript::Data r=escript::Data(res, what);
3065            return r;
3066        }
3067        return res;
3068    }
3069    
3070  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
3071  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 3104  inset=2*radius+1
3104  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
3105  that ripley has.  that ripley has.
3106  */  */
3107  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
3108  {  {
3109      if (m_numDim!=3)      unsigned int radius=0;  // these are only used by gaussian
3110        double sigma=0.5;
3111        
3112        unsigned int numvals=escript::DataTypes::noValues(shape);
3113        
3114        if (len(filter)==0)
3115      {      {
3116          throw RipleyException("Brick must be 3D.");      // nothing special required here yet
3117      }      }
3118      if (len(filter)!=3) {      else if (len(filter)==3)
         throw RipleyException("Unsupported random filter");  
     }  
     boost::python::extract<string> ex(filter[0]);  
     if (!ex.check() || (ex()!="gaussian"))  
3119      {      {
3120          throw RipleyException("Unsupported random filter");          boost::python::extract<string> ex(filter[0]);
3121            if (!ex.check() || (ex()!="gaussian"))
3122            {
3123                throw RipleyException("Unsupported random filter for Brick.");
3124            }
3125            boost::python::extract<unsigned int> ex1(filter[1]);
3126            if (!ex1.check())
3127            {
3128                throw RipleyException("Radius of gaussian filter must be a positive integer.");
3129            }
3130            radius=ex1();
3131            sigma=0.5;
3132            boost::python::extract<double> ex2(filter[2]);
3133            if (!ex2.check() || (sigma=ex2())<=0)
3134            {
3135                throw RipleyException("Sigma must be a positive floating point number.");
3136            }            
3137      }      }
3138      boost::python::extract<unsigned int> ex1(filter[1]);      else
     if (!ex1.check())  
3139      {      {
3140          throw RipleyException("Radius of gaussian filter must be a positive integer.");          throw RipleyException("Unsupported random filter");
3141      }      }
3142      unsigned int radius=ex1();  
3143      double sigma=0.5;      // number of points in the internal region
3144      boost::python::extract<double> ex2(filter[2]);      // that is, the ones we need smoothed versions of
3145      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  
3146      size_t ext[3];      size_t ext[3];
3147      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
3148      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3149      ext[2]=internal[2]+2*radius;    // for smoothing      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3150            
3151      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3152      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3153    
3154      if ((2*radius>=internal[0]) || (2*radius>=internal[1]) || (2*radius>=internal[2]))      if (2*radius>=internal[0]-4)
3155      {      {
3156          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");
3157      }      }
3158            if (2*radius>=internal[1]-4)
3159        {
3160            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3161        }
3162        if (2*radius>=internal[2]-4)
3163        {
3164            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3165        }    
3166        
3167      double* src=new double[ext[0]*ext[1]*ext[2]];      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3168      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3169            
     
3170  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3171          
3172        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3173        {
3174        // since the dimensions are equal for all ranks, this exception
3175        // will be thrown on all ranks
3176        throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3177    
3178        }
3179      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3180      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];
3181      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3182    #endif    
3183    
3184    /*    
3185        // if we wanted to test a repeating pattern
3186        size_t basex=0;
3187        size_t basey=0;
3188        size_t basez=0;
3189    #ifdef ESYS_MPI    
3190        basex=X*m_gNE[0]/m_NX[0];
3191        basey=Y*m_gNE[1]/m_NX[1];
3192        basez=Z*m_gNE[2]/m_NX[2];
3193            
3194    cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
3195        
3196    #endif    
3197        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3198    */
3199    
3200    #ifdef ESYS_MPI
3201    
3202    
3203    
3204      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);
3205      size_t inset=2*radius+1;          size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3206            // there is an overlap of two points both of which need to have "radius" points on either side.
3207            
3208      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
3209      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;  
3210      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3211            
3212      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);    
3213            
3214      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
3215      MPI_Status stats[50];      MPI_Status stats[50];
3216      short rused=0;      short rused=0;
3217            
# Line 3019  escript::Data Brick::randomFill(long see Line 3222  escript::Data Brick::randomFill(long see
3222      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3223            
3224            
3225      block.copyUsedFromBuffer(src);      block.copyAllToBuffer(src);
       
3226            
3227      int comserr=0;          int comserr=0;    
3228      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
3229      {      {
3230      message& m=incoms[i];          message& m=incoms[i];
3231      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++));
3232      block.setUsed(m.destbuffid);          block.setUsed(m.destbuffid);
3233      }      }
3234    
3235      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3236      {      {
3237      message& m=incoms[i];          message& m=outcoms[i];
3238      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++));
3239      }          }    
3240            
3241      if (!comserr)      if (!comserr)
# Line 3043  escript::Data Brick::randomFill(long see Line 3245  escript::Data Brick::randomFill(long see
3245    
3246      if (comserr)      if (comserr)
3247      {      {
3248      // 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.
3249      // 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.
3250      // 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
3251          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
3252      }      }
3253            
3254      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3255        
3256  #endif      #endif    
3257      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      
3258      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
     // don't need to check for exwrite because we just made it  
     escript::DataVector& dv=resdat.getExpandedVectorReference();  
     double* convolution=get3DGauss(radius, sigma);  
     for (size_t z=0;z<(internal[2]);++z)  
3259      {      {
3260      for (size_t y=0;y<(internal[1]);++y)            
3261      {          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3262          for (size_t x=0;x<(internal[0]);++x)          escript::Data resdat(0, shape, fs , true);
3263          {              // don't need to check for exwrite because we just made it
3264          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();
3265                
3266          }      
3267      }          // now we need to copy values over
3268      }          for (size_t z=0;z<(internal[2]);++z)
3269      delete[] convolution;          {
3270      delete[] src;              for (size_t y=0;y<(internal[1]);++y)    
3271      return resdat;              {
3272                    for (size_t x=0;x<(internal[0]);++x)
3273                    {
3274                        for (unsigned int i=0;i<numvals;++i)
3275                        {
3276                            dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3277                        }
3278                    }
3279                }
3280            }  
3281            delete[] src;
3282            return resdat;      
3283        }
3284        else        // filter enabled  
3285        {
3286        
3287            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3288            escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3289            // don't need to check for exwrite because we just made it
3290            escript::DataVector& dv=resdat.getExpandedVectorReference();
3291            double* convolution=get3DGauss(radius, sigma);
3292    
3293            for (size_t z=0;z<(internal[2]);++z)
3294            {
3295                for (size_t y=0;y<(internal[1]);++y)    
3296                {
3297                    for (size_t x=0;x<(internal[0]);++x)
3298                    {    
3299                        dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3300                
3301                    }
3302                }
3303            }
3304        
3305            delete[] convolution;
3306            delete[] src;
3307            return resdat;
3308        
3309        }
3310  }  }
3311    
3312    
# Line 3081  escript::Data Brick::randomFill(long see Line 3317  escript::Data Brick::randomFill(long see
3317  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const {
3318      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3319      //is the found element even owned by this rank      //is the found element even owned by this rank
3320        // (inside owned or shared elements but will map to an owned element)
3321      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3322          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3323                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
3324            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3325                    + m_dx[dim]/2.;
3326            if (min > coords[dim] || max < coords[dim]) {
3327              return NOT_MINE;              return NOT_MINE;
3328          }          }
3329      }      }
# Line 3102  int Brick::findNode(const double *coords Line 3342  int Brick::findNode(const double *coords
3342          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
3343      }      }
3344      //find the closest node      //find the closest node
3345      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
3346          double xdist = x - (ex + dx)*m_dx[0];          double xdist = x - (ex + dx)*m_dx[0];
3347          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
3348              double ydist = y - (ey + dy)*m_dx[1];              double ydist = y - (ey + dy)*m_dx[1];
3349              for (int dz = 0; dz < 2; dz++) {              for (int dz = 0; dz < 1; dz++) {
3350                  double zdist = z - (ez + dz)*m_dx[2];                  double zdist = z - (ez + dz)*m_dx[2];
3351                  double total = xdist*xdist + ydist*ydist + zdist*zdist;                  double total = xdist*xdist + ydist*ydist + zdist*zdist;
3352                  if (total < minDist) {                  if (total < minDist) {
3353                      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],
3354                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3355                      minDist = total;                      minDist = total;
3356                  }                  }
3357              }              }
3358          }          }
3359      }      }
3360        if (closest == NOT_MINE) {
3361            throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");
3362        }
3363      return closest;      return closest;
3364  }  }
3365    
3366  void Brick::setAssembler(std::string type, std::map<std::string,  void Brick::setAssembler(std::string type, std::map<std::string,
3367          escript::Data> constants) {          escript::Data> constants) {
3368      if (type.compare("WaveAssembler") == 0) {      if (type.compare("WaveAssembler") == 0) {
3369            if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3370                throw RipleyException("Domain already using a different custom assembler");
3371            assembler_type = WAVE_ASSEMBLER;
3372          delete assembler;          delete assembler;
3373          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);
3374        } else if (type.compare("LameAssembler") == 0) {
3375            if (assembler_type != LAME_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3376                throw RipleyException("Domain already using a different custom assembler");
3377            assembler_type = LAME_ASSEMBLER;
3378            delete assembler;
3379            assembler = new LameAssembler3D(this, m_dx, m_NX, m_NE, m_NN);
3380      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3381          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Brick does not support the"
3382                                  " requested assembler");                                  " requested assembler");
3383      }      }
3384  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26