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

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

  ViewVC Help
Powered by ViewVC 1.1.26