/[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 4860 by sshaw, Thu Apr 10 04:08:42 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,
63               const simap_t& tagnamestonums) :               const simap_t& tagnamestonums) :
64      RipleyDomain(3)      RipleyDomain(3)
65  {  {
66        if (n0 <= 0 || n1 <= 0 || n2 <= 0)
67            throw RipleyException("Number of elements in each spatial dimension "
68                    "must be positive");
69    
70      // ignore subdivision parameters for serial run      // ignore subdivision parameters for serial run
71      if (m_mpiInfo->size == 1) {      if (m_mpiInfo->size == 1) {
72          d0=1;          d0=1;
# Line 55  Brick::Brick(int n0, int n1, int n2, dou Line 74  Brick::Brick(int n0, int n1, int n2, dou
74          d2=1;          d2=1;
75      }      }
76      bool warn=false;      bool warn=false;
77      // if number of subdivisions is non-positive, try to subdivide by the same  
78      // ratio as the number of elements      std::vector<int> factors;
79      if (d0<=0 && d1<=0 && d2<=0) {      int ranks = m_mpiInfo->size;
80          warn=true;      int epr[3] = {n0,n1,n2};
81          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));      int d[3] = {d0,d1,d2};
82          d0=max(1, d0);      if (d0<=0 || d1<=0 || d2<=0) {
83          d1=max(1, (int)(d0*n1/(float)n0));          for (int i = 0; i < 3; i++) {
84          d2=m_mpiInfo->size/(d0*d1);              if (d[i] < 1) {
85          if (d0*d1*d2 != m_mpiInfo->size) {                  d[i] = 1;
86              // ratios not the same so leave "smallest" side undivided and try                  continue;
87              // dividing 2 sides only              }
88              if (n0>=n1) {              epr[i] = -1; // can no longer be max
89                  if (n1>=n2) {              if (ranks % d[i] != 0) {
90                      d0=d1=0;                  throw RipleyException("Invalid number of spatial subdivisions");
91                      d2=1;              }
92                  } else {              //remove
93                      d0=d2=0;              ranks /= d[i];
94                      d1=1;          }
95                  }          factorise(factors, ranks);
96              } else {          if (factors.size() != 0) {
97                  if (n0>=n2) {              warn = true;
98                      d0=d1=0;          }
99                      d2=1;      }
100                  } else {      while (factors.size() > 0) {
101                      d0=1;          int i = indexOfMax(epr[0],epr[1],epr[2]);
102                      d1=d2=0;          int f = factors.back();
103                  }          factors.pop_back();
104              }          d[i] *= f;
105          }          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);  
106      }      }
107        d0 = d[0]; d1 = d[1]; d2 = d[2];
108    
109      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
110      // among number of ranks      // among number of ranks
111      if (d0*d1*d2 != m_mpiInfo->size)      if (d0*d1*d2 != m_mpiInfo->size){
112          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
113        }
114      if (warn) {      if (warn) {
115          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
116              << 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 207  Brick::Brick(int n0, int n1, int n2, dou
207    
208  Brick::~Brick()  Brick::~Brick()
209  {  {
     Paso_SystemMatrixPattern_free(m_pattern);  
     Paso_Connector_free(m_connector);  
210      delete assembler;      delete assembler;
211  }  }
212    
# Line 411  void Brick::readNcGrid(escript::Data& ou Line 376  void Brick::readNcGrid(escript::Data& ou
376  #endif  #endif
377  }  }
378    
379    #ifdef USE_BOOSTIO
380    void Brick::readBinaryGridFromZipped(escript::Data& out, string filename,
381                               const ReaderParameters& params) const
382    {
383        // the mapping is not universally correct but should work on our
384        // supported platforms
385        switch (params.dataType) {
386            case DATATYPE_INT32:
387                readBinaryGridZippedImpl<int>(out, filename, params);
388                break;
389            case DATATYPE_FLOAT32:
390                readBinaryGridZippedImpl<float>(out, filename, params);
391                break;
392            case DATATYPE_FLOAT64:
393                readBinaryGridZippedImpl<double>(out, filename, params);
394                break;
395            default:
396                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
397        }
398    }
399    #endif
400    
401  void Brick::readBinaryGrid(escript::Data& out, string filename,  void Brick::readBinaryGrid(escript::Data& out, string filename,
402                             const ReaderParameters& params) const                             const ReaderParameters& params) const
403  {  {
# Line 547  void Brick::readBinaryGridImpl(escript:: Line 534  void Brick::readBinaryGridImpl(escript::
534      f.close();      f.close();
535  }  }
536    
537    #ifdef USE_BOOSTIO
538    template<typename ValueType>
539    void Brick::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
540                                   const ReaderParameters& params) const
541    {
542        // check destination function space
543        int myN0, myN1, myN2;
544        if (out.getFunctionSpace().getTypeCode() == Nodes) {
545            myN0 = m_NN[0];
546            myN1 = m_NN[1];
547            myN2 = m_NN[2];
548        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
549                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
550            myN0 = m_NE[0];
551            myN1 = m_NE[1];
552            myN2 = m_NE[2];
553        } else
554            throw RipleyException("readBinaryGridFromZipped(): invalid function space for output data object");
555    
556        if (params.first.size() != 3)
557            throw RipleyException("readBinaryGridFromZipped(): argument 'first' must have 3 entries");
558    
559        if (params.numValues.size() != 3)
560            throw RipleyException("readBinaryGridFromZipped(): argument 'numValues' must have 3 entries");
561    
562        if (params.multiplier.size() != 3)
563            throw RipleyException("readBinaryGridFromZipped(): argument 'multiplier' must have 3 entries");
564        for (size_t i=0; i<params.multiplier.size(); i++)
565            if (params.multiplier[i]<1)
566                throw RipleyException("readBinaryGridFromZipped(): all multipliers must be positive");
567    
568        // check file existence and size
569        ifstream f(filename.c_str(), ifstream::binary);
570        if (f.fail()) {
571            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
572        }
573        f.seekg(0, ios::end);
574        const int numComp = out.getDataPointSize();
575        int filesize = f.tellg();
576        f.seekg(0, ios::beg);
577        std::vector<char> compressed(filesize);
578        f.read((char*)&compressed[0], filesize);
579        f.close();
580        std::vector<char> decompressed = unzip(compressed);
581        filesize = decompressed.size();
582        const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
583        if (filesize < reqsize) {
584            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
585        }
586    
587        // check if this rank contributes anything
588        if (params.first[0] >= m_offset[0]+myN0 ||
589                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
590                params.first[1] >= m_offset[1]+myN1 ||
591                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
592                params.first[2] >= m_offset[2]+myN2 ||
593                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
594            return;
595        }
596    
597        // now determine how much this rank has to write
598    
599        // first coordinates in data object to write to
600        const int first0 = max(0, params.first[0]-m_offset[0]);
601        const int first1 = max(0, params.first[1]-m_offset[1]);
602        const int first2 = max(0, params.first[2]-m_offset[2]);
603        // indices to first value in file
604        const int idx0 = max(0, m_offset[0]-params.first[0]);
605        const int idx1 = max(0, m_offset[1]-params.first[1]);
606        const int idx2 = max(0, m_offset[2]-params.first[2]);
607        // number of values to read
608        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
609        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
610        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
611    
612        out.requireWrite();
613        vector<ValueType> values(num0*numComp);
614        const int dpp = out.getNumDataPointsPerSample();
615    
616        for (int z=0; z<num2; z++) {
617            for (int y=0; y<num1; y++) {
618                const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
619                                 +(idx2+z)*params.numValues[0]*params.numValues[1]);
620                memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
621                
622                for (int x=0; x<num0; x++) {
623                    const int baseIndex = first0+x*params.multiplier[0]
624                                         +(first1+y*params.multiplier[1])*myN0
625                                         +(first2+z*params.multiplier[2])*myN0*myN1;
626                    for (int m2=0; m2<params.multiplier[2]; m2++) {
627                        for (int m1=0; m1<params.multiplier[1]; m1++) {
628                            for (int m0=0; m0<params.multiplier[0]; m0++) {
629                                const int dataIndex = baseIndex+m0
630                                               +m1*myN0
631                                               +m2*myN0*myN1;
632                                double* dest = out.getSampleDataRW(dataIndex);
633                                for (int c=0; c<numComp; c++) {
634                                    ValueType val = values[x*numComp+c];
635    
636                                    if (params.byteOrder != BYTEORDER_NATIVE) {
637                                        char* cval = reinterpret_cast<char*>(&val);
638                                        // this will alter val!!
639                                        byte_swap32(cval);
640                                    }
641                                    if (!std::isnan(val)) {
642                                        for (int q=0; q<dpp; q++) {
643                                            *dest++ = static_cast<double>(val);
644                                        }
645                                    }
646                                }
647                            }
648                        }
649                    }
650                }
651            }
652        }
653    }
654    #endif
655    
656  void Brick::writeBinaryGrid(const escript::Data& in, string filename,  void Brick::writeBinaryGrid(const escript::Data& in, string filename,
657                              int byteOrder, int dataType) const                              int byteOrder, int dataType) const
658  {  {
# Line 786  const int* Brick::borrowSampleReferenceI Line 892  const int* Brick::borrowSampleReferenceI
892          case FaceElements:          case FaceElements:
893          case ReducedFaceElements:          case ReducedFaceElements:
894              return &m_faceId[0];              return &m_faceId[0];
895            case Points:
896                return &m_diracPointNodeIDs[0];
897          default:          default:
898              break;              break;
899      }      }
# Line 1944  void Brick::nodesToDOF(escript::Data& ou Line 2052  void Brick::nodesToDOF(escript::Data& ou
2052  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const
2053  {  {
2054      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
2055      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2056      // 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
2057      const_cast<escript::Data*>(&in)->expand();      const_cast<escript::Data*>(&in)->expand();
2058      Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));      coupler->startCollect(in.getSampleDataRO(0));
2059    
2060      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2061      out.requireWrite();      out.requireWrite();
2062      const double* buffer = Paso_Coupler_finishCollect(coupler);      const double* buffer = coupler->finishCollect();
2063    
2064  #pragma omp parallel for  #pragma omp parallel for
2065      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 2068  void Brick::dofToNodes(escript::Data& ou
2068                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
2069          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
2070      }      }
     Paso_Coupler_free(coupler);  
2071  }  }
2072    
2073  //private  //private
# Line 1981  void Brick::populateSampleIds() Line 2088  void Brick::populateSampleIds()
2088          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2089      }      }
2090      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2091      m_nodeId.resize(getNumNodes());      
2092      m_dofId.resize(numDOF);      try {
2093      m_elementId.resize(getNumElements());          m_nodeId.resize(getNumNodes());
2094            m_dofId.resize(numDOF);
2095            m_elementId.resize(getNumElements());
2096        } catch (const std::length_error& le) {
2097            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2098        }
2099        
2100      // populate face element counts      // populate face element counts
2101      //left      //left
2102      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2228  void Brick::createPattern() Line 2340  void Brick::createPattern()
2340      RankVector neighbour;      RankVector neighbour;
2341      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
2342      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
2343      int numShared=0;      int numShared=0, expectedShared=0;;
2344      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
2345      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];
2346      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 2354  void Brick::createPattern()
2354                  const int nx=x+i0;                  const int nx=x+i0;
2355                  const int ny=y+i1;                  const int ny=y+i1;
2356                  const int nz=z+i2;                  const int nz=z+i2;
2357                    if (!(nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2])) {
2358                        continue;
2359                    }
2360                    if (i0==0 && i1==0)
2361                        expectedShared += nDOF0*nDOF1;
2362                    else if (i0==0 && i2==0)
2363                        expectedShared += nDOF0*nDOF2;
2364                    else if (i1==0 && i2==0)
2365                        expectedShared += nDOF1*nDOF2;
2366                    else if (i0==0)
2367                        expectedShared += nDOF0;
2368                    else if (i1==0)
2369                        expectedShared += nDOF1;
2370                    else if (i2==0)
2371                        expectedShared += nDOF2;
2372                    else
2373                        expectedShared++;
2374                }
2375            }
2376        }
2377        
2378        vector<IndexVector> rowIndices(expectedShared);
2379        
2380        for (int i2=-1; i2<2; i2++) {
2381            for (int i1=-1; i1<2; i1++) {
2382                for (int i0=-1; i0<2; i0++) {
2383                    // skip this rank
2384                    if (i0==0 && i1==0 && i2==0)
2385                        continue;
2386                    // location of neighbour rank
2387                    const int nx=x+i0;
2388                    const int ny=y+i1;
2389                    const int nz=z+i2;
2390                  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]) {
2391                      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);
2392                      if (i0==0 && i1==0) {                      if (i0==0 && i1==0) {
# Line 2257  void Brick::createPattern() Line 2402  void Brick::createPattern()
2402                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2403                                  if (j>0) {                                  if (j>0) {
2404                                      if (i>0)                                      if (i>0)
2405                                          colIndices[firstDOF+j-1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0, numShared);
2406                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2407                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2408                                          colIndices[firstDOF+j-1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0, numShared);
2409                                  }                                  }
2410                                  if (i>0)                                  if (i>0)
2411                                      colIndices[firstDOF+j-nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0, numShared);
2412                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2413                                  if (i<nDOF1-1)                                  if (i<nDOF1-1)
2414                                      colIndices[firstDOF+j+nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0, numShared);
2415                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2416                                      if (i>0)                                      if (i>0)
2417                                          colIndices[firstDOF+j+1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0, numShared);
2418                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2419                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2420                                          colIndices[firstDOF+j+1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0, numShared);
2421                                  }                                  }
2422                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2423                              }                              }
# Line 2291  void Brick::createPattern() Line 2436  void Brick::createPattern()
2436                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2437                                  if (j>0) {                                  if (j>0) {
2438                                      if (i>0)                                      if (i>0)
2439                                          colIndices[firstDOF+j-1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0*nDOF1, numShared);
2440                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2441                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2442                                          colIndices[firstDOF+j-1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0*nDOF1, numShared);
2443                                  }                                  }
2444                                  if (i>0)                                  if (i>0)
2445                                      colIndices[firstDOF+j-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0*nDOF1, numShared);
2446                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2447                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2448                                      colIndices[firstDOF+j+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0*nDOF1, numShared);
2449                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2450                                      if (i>0)                                      if (i>0)
2451                                          colIndices[firstDOF+j+1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0*nDOF1, numShared);
2452                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2453                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2454                                          colIndices[firstDOF+j+1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0*nDOF1, numShared);
2455                                  }                                  }
2456                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2457                              }                              }
# Line 2325  void Brick::createPattern() Line 2470  void Brick::createPattern()
2470                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2471                                  if (j>0) {                                  if (j>0) {
2472                                      if (i>0)                                      if (i>0)
2473                                          colIndices[firstDOF+(j-1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0-nDOF0*nDOF1, numShared);
2474                                      colIndices[firstDOF+(j-1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0, numShared);
2475                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2476                                          colIndices[firstDOF+(j-1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0+nDOF0*nDOF1, numShared);
2477                                  }                                  }
2478                                  if (i>0)                                  if (i>0)
2479                                      colIndices[firstDOF+j*nDOF0-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0-nDOF0*nDOF1, numShared);
2480                                  colIndices[firstDOF+j*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0, numShared);
2481                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2482                                      colIndices[firstDOF+j*nDOF0+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0+nDOF0*nDOF1, numShared);
2483                                  if (j<nDOF1-1) {                                  if (j<nDOF1-1) {
2484                                      if (i>0)                                      if (i>0)
2485                                          colIndices[firstDOF+(j+1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0-nDOF0*nDOF1, numShared);
2486                                      colIndices[firstDOF+(j+1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0, numShared);
2487                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2488                                          colIndices[firstDOF+(j+1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0+nDOF0*nDOF1, numShared);
2489                                  }                                  }
2490                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2491                              }                              }
# Line 2356  void Brick::createPattern() Line 2501  void Brick::createPattern()
2501                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
2502                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2503                              if (i>0)                              if (i>0)
2504                                  colIndices[firstDOF+i-1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
2505                              colIndices[firstDOF+i].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
2506                              if (i<nDOF0-1)                              if (i<nDOF0-1)
2507                                  colIndices[firstDOF+i+1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
2508                              m_dofMap[firstNode+i]=numDOF+numShared;                              m_dofMap[firstNode+i]=numDOF+numShared;
2509                          }                          }
2510                      } else if (i1==0) {                      } else if (i1==0) {
# Line 2374  void Brick::createPattern() Line 2519  void Brick::createPattern()
2519                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
2520                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2521                              if (i>0)                              if (i>0)
2522                                  colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
2523                              colIndices[firstDOF+i*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
2524                              if (i<nDOF1-1)                              if (i<nDOF1-1)
2525                                  colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
2526                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2527                          }                          }
2528                      } else if (i2==0) {                      } else if (i2==0) {
# Line 2392  void Brick::createPattern() Line 2537  void Brick::createPattern()
2537                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2538                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2539                              if (i>0)                              if (i>0)
2540                                  colIndices[firstDOF+(i-1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0*nDOF1, numShared);
2541                              colIndices[firstDOF+i*nDOF0*nDOF1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0*nDOF1, numShared);
2542                              if (i<nDOF2-1)                              if (i<nDOF2-1)
2543                                  colIndices[firstDOF+(i+1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0*nDOF1, numShared);
2544                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2545                          }                          }
2546                      } else {                      } else {
# Line 2409  void Brick::createPattern() Line 2554  void Brick::createPattern()
2554                          offsetInShared.push_back(offsetInShared.back()+1);                          offsetInShared.push_back(offsetInShared.back()+1);
2555                          sendShared.push_back(dof);                          sendShared.push_back(dof);
2556                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
2557                          colIndices[dof].push_back(numShared);                          doublyLink(colIndices, rowIndices, dof, numShared);
2558                          m_dofMap[node]=numDOF+numShared;                          m_dofMap[node]=numDOF+numShared;
2559                          ++numShared;                          ++numShared;
2560                      }                      }
# Line 2418  void Brick::createPattern() Line 2563  void Brick::createPattern()
2563          }          }
2564      }      }
2565    
2566    #pragma omp parallel for
2567        for (int i = 0; i < numShared; i++) {
2568            std::sort(rowIndices[i].begin(), rowIndices[i].end());
2569        }
2570    
2571      // create connector      // create connector
2572      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
2573              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2574              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2575      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2576              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2577              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2578      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);  
2579    
2580      // create main and couple blocks      // create main and couple blocks
2581      Paso_Pattern *mainPattern = createMainPattern();      paso::Pattern_ptr mainPattern = createMainPattern();
2582      Paso_Pattern *colPattern, *rowPattern;      paso::Pattern_ptr colPattern, rowPattern;
2583      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, colPattern, rowPattern);
2584    
2585      // allocate paso distribution      // allocate paso distribution
2586      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2587              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2588    
2589      // finally create the system matrix      // finally create the system matrix
2590      m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,      m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2591              distribution, distribution, mainPattern, colPattern, rowPattern,              distribution, distribution, mainPattern, colPattern, rowPattern,
2592              m_connector, m_connector);              m_connector, m_connector));
   
     Paso_Distribution_free(distribution);  
2593    
2594      // useful debug output      // useful debug output
2595      /*      /*
# Line 2502  void Brick::createPattern() Line 2648  void Brick::createPattern()
2648          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2649      }      }
2650      */      */
   
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colPattern);  
     Paso_Pattern_free(rowPattern);  
2651  }  }
2652    
2653  //private  //private
2654  void Brick::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,  void Brick::addToMatrixAndRHS(paso::SystemMatrix_ptr S, escript::Data& F,
2655           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2656           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2657  {  {
# Line 2858  void Brick::interpolateNodesOnFaces(escr Line 3000  void Brick::interpolateNodesOnFaces(escr
3000    
3001  namespace  namespace
3002  {  {
3003      // Calculates a guassian blur colvolution matrix for 3D      // Calculates a gaussian blur convolution matrix for 3D
3004      // See wiki article on the subject      // See wiki article on the subject
3005      double* get3DGauss(unsigned radius, double sigma)      double* get3DGauss(unsigned radius, double sigma)
3006      {      {
3007          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)];
3008          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);
3009      double total=0;          double total=0;
3010      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
3011      for (int z=-r;z<=r;++z)          for (int z=-r;z<=r;++z)
3012      {          {
3013          for (int y=-r;y<=r;++y)              for (int y=-r;y<=r;++y)
3014          {              {
3015          for (int x=-r;x<=r;++x)                  for (int x=-r;x<=r;++x)
3016          {                          {        
3017              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));
3018              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)];    
3019          }                  }
3020          }              }
3021      }          }
3022      double invtotal=1/total;          double invtotal=1/total;
3023      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)
3024      {          {
3025          arr[p]*=invtotal;              arr[p]*=invtotal;
3026      }          }
3027      return arr;          return arr;
3028      }      }
3029            
3030      // applies conv to source to get a point.      // applies conv to source to get a point.
# Line 2890  namespace Line 3032  namespace
3032      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)
3033      {      {
3034          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3035      size_t sbase=bx+by*width+bz*width*height;          size_t sbase=bx+by*width+bz*width*height;
3036      double result=0;          double result=0;
3037      for (int z=0;z<2*radius+1;++z)          for (int z=0;z<2*radius+1;++z)
3038      {          {
3039          for (int y=0;y<2*radius+1;++y)              for (int y=0;y<2*radius+1;++y)
3040          {                  {    
3041          for (int x=0;x<2*radius+1;++x)                  for (int x=0;x<2*radius+1;++x)
3042          {                  {
3043              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];
3044          }                  }
3045          }              }
3046      }          }
3047            // use this line for "pass-though" (return the centre point value)
3048    //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3049          return result;                return result;      
3050      }      }
3051  }  }
3052    
3053    /* This is a wrapper for filtered (and non-filtered) randoms
3054     * For detailed doco see randomFillWorker
3055    */
3056    escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
3057           const escript::FunctionSpace& what,
3058           long seed, const boost::python::tuple& filter) const
3059    {
3060        int numvals=escript::DataTypes::noValues(shape);
3061        if (len(filter)>0 && (numvals!=1))
3062        {
3063            throw RipleyException("Ripley only supports filters for scalar data.");
3064        }
3065        escript::Data res=randomFillWorker(shape, seed, filter);
3066        if (res.getFunctionSpace()!=what)
3067        {
3068            escript::Data r=escript::Data(res, what);
3069            return r;
3070        }
3071        return res;
3072    }
3073    
3074  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
3075  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 3108  inset=2*radius+1
3108  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
3109  that ripley has.  that ripley has.
3110  */  */
3111  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
3112  {  {
3113      if (m_numDim!=3)      unsigned int radius=0;  // these are only used by gaussian
3114        double sigma=0.5;
3115        
3116        unsigned int numvals=escript::DataTypes::noValues(shape);
3117        
3118        if (len(filter)==0)
3119      {      {
3120          throw RipleyException("Brick must be 3D.");      // nothing special required here yet
     }  
     if (len(filter)!=3) {  
         throw RipleyException("Unsupported random filter");  
3121      }      }
3122      boost::python::extract<string> ex(filter[0]);      else if (len(filter)==3)
     if (!ex.check() || (ex()!="gaussian"))  
3123      {      {
3124          throw RipleyException("Unsupported random filter");          boost::python::extract<string> ex(filter[0]);
3125            if (!ex.check() || (ex()!="gaussian"))
3126            {
3127                throw RipleyException("Unsupported random filter for Brick.");
3128            }
3129            boost::python::extract<unsigned int> ex1(filter[1]);
3130            if (!ex1.check())
3131            {
3132                throw RipleyException("Radius of gaussian filter must be a positive integer.");
3133            }
3134            radius=ex1();
3135            sigma=0.5;
3136            boost::python::extract<double> ex2(filter[2]);
3137            if (!ex2.check() || (sigma=ex2())<=0)
3138            {
3139                throw RipleyException("Sigma must be a positive floating point number.");
3140            }            
3141      }      }
3142      boost::python::extract<unsigned int> ex1(filter[1]);      else
     if (!ex1.check())  
3143      {      {
3144          throw RipleyException("Radius of gaussian filter must be a positive integer.");          throw RipleyException("Unsupported random filter");
3145      }      }
3146      unsigned int radius=ex1();  
3147      double sigma=0.5;      // number of points in the internal region
3148      boost::python::extract<double> ex2(filter[2]);      // that is, the ones we need smoothed versions of
3149      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  
3150      size_t ext[3];      size_t ext[3];
3151      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
3152      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3153      ext[2]=internal[2]+2*radius;    // for smoothing      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3154            
3155      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3156      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3157    
3158      if ((2*radius>=internal[0]) || (2*radius>=internal[1]) || (2*radius>=internal[2]))      if (2*radius>=internal[0]-4)
3159      {      {
3160          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");
3161      }      }
3162            if (2*radius>=internal[1]-4)
3163        {
3164            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3165        }
3166        if (2*radius>=internal[2]-4)
3167        {
3168            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3169        }    
3170        
3171      double* src=new double[ext[0]*ext[1]*ext[2]];      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3172      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3173            
     
3174  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3175          
3176        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3177        {
3178        // since the dimensions are equal for all ranks, this exception
3179        // will be thrown on all ranks
3180        throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3181    
3182        }
3183      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3184      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];
3185      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3186    #endif    
3187    
3188    /*    
3189        // if we wanted to test a repeating pattern
3190        size_t basex=0;
3191        size_t basey=0;
3192        size_t basez=0;
3193    #ifdef ESYS_MPI    
3194        basex=X*m_gNE[0]/m_NX[0];
3195        basey=Y*m_gNE[1]/m_NX[1];
3196        basez=Z*m_gNE[2]/m_NX[2];
3197            
3198    cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
3199        
3200    #endif    
3201        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3202    */
3203    
3204    #ifdef ESYS_MPI
3205    
3206    
3207    
3208      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);
3209      size_t inset=2*radius+1;          size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3210            // there is an overlap of two points both of which need to have "radius" points on either side.
3211            
3212      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
3213      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;  
3214      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3215            
3216      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);    
3217            
3218      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
3219      MPI_Status stats[50];      MPI_Status stats[50];
3220      short rused=0;      short rused=0;
3221            
# Line 3019  escript::Data Brick::randomFill(long see Line 3226  escript::Data Brick::randomFill(long see
3226      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3227            
3228            
3229      block.copyUsedFromBuffer(src);      block.copyAllToBuffer(src);
       
3230            
3231      int comserr=0;          int comserr=0;    
3232      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
3233      {      {
3234      message& m=incoms[i];          message& m=incoms[i];
3235      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++));
3236      block.setUsed(m.destbuffid);          block.setUsed(m.destbuffid);
3237      }      }
3238    
3239      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3240      {      {
3241      message& m=incoms[i];          message& m=outcoms[i];
3242      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++));
3243      }          }    
3244            
3245      if (!comserr)      if (!comserr)
# Line 3043  escript::Data Brick::randomFill(long see Line 3249  escript::Data Brick::randomFill(long see
3249    
3250      if (comserr)      if (comserr)
3251      {      {
3252      // 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.
3253      // 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.
3254      // 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
3255          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
3256      }      }
3257            
3258      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3259        
3260  #endif      #endif    
3261      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      
3262      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
3263      // don't need to check for exwrite because we just made it      {
3264      escript::DataVector& dv=resdat.getExpandedVectorReference();        
3265      double* convolution=get3DGauss(radius, sigma);          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3266      for (size_t z=0;z<(internal[2]);++z)          escript::Data resdat(0, shape, fs , true);
3267            // don't need to check for exwrite because we just made it
3268            escript::DataVector& dv=resdat.getExpandedVectorReference();
3269        
3270        
3271            // now we need to copy values over
3272            for (size_t z=0;z<(internal[2]);++z)
3273            {
3274                for (size_t y=0;y<(internal[1]);++y)    
3275                {
3276                    for (size_t x=0;x<(internal[0]);++x)
3277                    {
3278                        for (unsigned int i=0;i<numvals;++i)
3279                        {
3280                            dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3281                        }
3282                    }
3283                }
3284            }  
3285            delete[] src;
3286            return resdat;      
3287        }
3288        else        // filter enabled  
3289      {      {
3290      for (size_t y=0;y<(internal[1]);++y)          
3291      {          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3292          for (size_t x=0;x<(internal[0]);++x)          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3293          {              // don't need to check for exwrite because we just made it
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]);          escript::DataVector& dv=resdat.getExpandedVectorReference();
3295                    double* convolution=get3DGauss(radius, sigma);
3296          }  
3297      }          for (size_t z=0;z<(internal[2]);++z)
3298      }          {
3299      delete[] convolution;              for (size_t y=0;y<(internal[1]);++y)    
3300      delete[] src;              {
3301      return resdat;                  for (size_t x=0;x<(internal[0]);++x)
3302                    {    
3303                        dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3304                
3305                    }
3306                }
3307            }
3308        
3309            delete[] convolution;
3310            delete[] src;
3311            return resdat;
3312        
3313        }
3314  }  }
3315    
3316    
# Line 3081  escript::Data Brick::randomFill(long see Line 3321  escript::Data Brick::randomFill(long see
3321  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const {
3322      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3323      //is the found element even owned by this rank      //is the found element even owned by this rank
3324        // (inside owned or shared elements but will map to an owned element)
3325      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3326          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3327                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
3328            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3329                    + m_dx[dim]/2.;
3330            if (min > coords[dim] || max < coords[dim]) {
3331              return NOT_MINE;              return NOT_MINE;
3332          }          }
3333      }      }
# Line 3091  int Brick::findNode(const double *coords Line 3335  int Brick::findNode(const double *coords
3335      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3336      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3337      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3338        
3339        //check if the point is even inside the domain
3340        if (x < 0 || y < 0 || z < 0
3341                || x > m_length[0] || y > m_length[1] || z > m_length[2])
3342            return NOT_MINE;
3343            
3344      // distance in elements      // distance in elements
3345      int ex = (int) floor(x / m_dx[0]);      int ex = (int) floor(x / m_dx[0]);
3346      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 3352  int Brick::findNode(const double *coords
3352          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
3353      }      }
3354      //find the closest node      //find the closest node
3355      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
3356          double xdist = x - (ex + dx)*m_dx[0];          double xdist = x - (ex + dx)*m_dx[0];
3357          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
3358              double ydist = y - (ey + dy)*m_dx[1];              double ydist = y - (ey + dy)*m_dx[1];
3359              for (int dz = 0; dz < 2; dz++) {              for (int dz = 0; dz < 1; dz++) {
3360                  double zdist = z - (ez + dz)*m_dx[2];                  double zdist = z - (ez + dz)*m_dx[2];
3361                  double total = xdist*xdist + ydist*ydist + zdist*zdist;                  double total = xdist*xdist + ydist*ydist + zdist*zdist;
3362                  if (total < minDist) {                  if (total < minDist) {
3363                      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],
3364                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3365                      minDist = total;                      minDist = total;
3366                  }                  }
3367              }              }
3368          }          }
3369      }      }
3370        if (closest == NOT_MINE) {
3371            throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");
3372        }
3373      return closest;      return closest;
3374  }  }
3375    
3376  void Brick::setAssembler(std::string type, std::map<std::string,  void Brick::setAssembler(std::string type, std::map<std::string,
3377          escript::Data> constants) {          escript::Data> constants) {
3378      if (type.compare("WaveAssembler") == 0) {      if (type.compare("WaveAssembler") == 0) {
3379            if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3380                throw RipleyException("Domain already using a different custom assembler");
3381            assembler_type = WAVE_ASSEMBLER;
3382          delete assembler;          delete assembler;
3383          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);
3384        } else if (type.compare("LameAssembler") == 0) {
3385            if (assembler_type != LAME_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3386                throw RipleyException("Domain already using a different custom assembler");
3387            assembler_type = LAME_ASSEMBLER;
3388            delete assembler;
3389            assembler = new LameAssembler3D(this, m_dx, m_NX, m_NE, m_NN);
3390      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3391          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Brick does not support the"
3392                                  " requested assembler");                                  " requested assembler");
3393      }      }
3394  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26