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

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

  ViewVC Help
Powered by ViewVC 1.1.26