/[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 4638 by jfenwick, Thu Jan 30 06:25:10 2014 UTC revision 4816 by caltinay, Fri Mar 28 06:16:02 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 33  Line 36 
36    
37  #include <iomanip>  #include <iomanip>
38    
39    #include "esysUtils/EsysRandom.h"
40    #include "blocktools.h"
41    
42    
43  using namespace std;  using namespace std;
44  using esysUtils::FileWriter;  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 51  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 236  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 407  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 543  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 782  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 1940  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 1956  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 2224  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 2238  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 2253  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 2287  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 2321  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 2352  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 2370  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 2388  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 2405  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 2414  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_ptr snd_shcomp(new paso::SharedComponents(
2567              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2568              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2569      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
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);
     Paso_SharedComponents_free(snd_shcomp);  
     Paso_SharedComponents_free(rcv_shcomp);  
2573    
2574      // create main and couple blocks      // create main and couple blocks
2575      Paso_Pattern *mainPattern = createMainPattern();      paso::Pattern *mainPattern = createMainPattern();
2576      Paso_Pattern *colPattern, *rowPattern;      paso::Pattern *colPattern, *rowPattern;
2577      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, &colPattern, &rowPattern);
2578    
2579      // allocate paso distribution      // allocate paso distribution
2580      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2581              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2582    
2583      // finally create the system matrix      // finally create the system matrix
2584      m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,      m_pattern = new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2585              distribution, distribution, mainPattern, colPattern, rowPattern,              distribution, distribution, mainPattern, colPattern, rowPattern,
2586              m_connector, m_connector);              m_connector, m_connector);
2587    
     Paso_Distribution_free(distribution);  
   
2588      // useful debug output      // useful debug output
2589      /*      /*
2590      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
# Line 2499  void Brick::createPattern() Line 2643  void Brick::createPattern()
2643      }      }
2644      */      */
2645    
2646      Paso_Pattern_free(mainPattern);      paso::Pattern_free(mainPattern);
2647      Paso_Pattern_free(colPattern);      paso::Pattern_free(colPattern);
2648      Paso_Pattern_free(rowPattern);      paso::Pattern_free(rowPattern);
2649  }  }
2650    
2651  //private  //private
# Line 2852  void Brick::interpolateNodesOnFaces(escr Line 2996  void Brick::interpolateNodesOnFaces(escr
2996      }      }
2997  }  }
2998    
2999    namespace
3000    {
3001        // Calculates a gaussian blur convolution matrix for 3D
3002        // See wiki article on the subject
3003        double* get3DGauss(unsigned radius, double sigma)
3004        {
3005            double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];
3006            double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);
3007            double total=0;
3008            int r=static_cast<int>(radius);
3009            for (int z=-r;z<=r;++z)
3010            {
3011                for (int y=-r;y<=r;++y)
3012                {
3013                    for (int x=-r;x<=r;++x)
3014                    {        
3015                        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));
3016                        total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];    
3017                    }
3018                }
3019            }
3020            double invtotal=1/total;
3021            for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
3022            {
3023                arr[p]*=invtotal;
3024            }
3025            return arr;
3026        }
3027        
3028        // applies conv to source to get a point.
3029        // (xp, yp) are the coords in the source matrix not the destination matrix
3030        double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)
3031        {
3032            size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3033            size_t sbase=bx+by*width+bz*width*height;
3034            double result=0;
3035            for (int z=0;z<2*radius+1;++z)
3036            {
3037                for (int y=0;y<2*radius+1;++y)
3038                {    
3039                    for (int x=0;x<2*radius+1;++x)
3040                    {
3041                        result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];
3042                    }
3043                }
3044            }
3045            // use this line for "pass-though" (return the centre point value)
3046    //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3047            return result;      
3048        }
3049    }
3050    
3051    /* This is a wrapper for filtered (and non-filtered) randoms
3052     * For detailed doco see randomFillWorker
3053    */
3054    escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
3055           const escript::FunctionSpace& what,
3056           long seed, const boost::python::tuple& filter) const
3057    {
3058        int numvals=escript::DataTypes::noValues(shape);
3059        if (len(filter)>0 && (numvals!=1))
3060        {
3061            throw RipleyException("Ripley only supports filters for scalar data.");
3062        }
3063        escript::Data res=randomFillWorker(shape, seed, filter);
3064        if (res.getFunctionSpace()!=what)
3065        {
3066            escript::Data r=escript::Data(res, what);
3067            return r;
3068        }
3069        return res;
3070    }
3071    
3072    /* This routine produces a Data object filled with smoothed random data.
3073    The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.
3074    A parameter radius  gives the size of the stencil used for the smoothing.
3075    A point on the left hand edge for example, will still require `radius` extra points to the left
3076    in order to complete the stencil.
3077    
3078    All local calculation is done on an array called `src`, with
3079    dimensions = ext[0] * ext[1] *ext[2].
3080    Where ext[i]= internal[i]+2*radius.
3081    
3082    Now for MPI there is overlap to deal with. We need to share both the overlapping
3083    values themselves but also the external region.
3084    
3085    In a hypothetical 1-D case:
3086    
3087    
3088    1234567
3089    would be split into two ranks thus:
3090    123(4)  (4)567     [4 being a shared element]
3091    
3092    If the radius is 2.   There will be padding elements on the outside:
3093    
3094    pp123(4)  (4)567pp
3095    
3096    To ensure that 4 can be correctly computed on both ranks, values from the other rank
3097    need to be known.
3098    
3099    pp123(4)56   23(4)567pp
3100    
3101    Now in our case, we wout set all the values 23456 on the left rank and send them to the
3102    right hand rank.
3103    
3104    So the edges _may_ need to be shared at a distance `inset` from all boundaries.
3105    inset=2*radius+1    
3106    This is to ensure that values at distance `radius` from the shared/overlapped element
3107    that ripley has.
3108    */
3109    escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
3110    {
3111        unsigned int radius=0;  // these are only used by gaussian
3112        double sigma=0.5;
3113        
3114        unsigned int numvals=escript::DataTypes::noValues(shape);
3115        
3116        if (len(filter)==0)
3117        {
3118        // nothing special required here yet
3119        }
3120        else if (len(filter)==3)
3121        {
3122            boost::python::extract<string> ex(filter[0]);
3123            if (!ex.check() || (ex()!="gaussian"))
3124            {
3125                throw RipleyException("Unsupported random filter for Brick.");
3126            }
3127            boost::python::extract<unsigned int> ex1(filter[1]);
3128            if (!ex1.check())
3129            {
3130                throw RipleyException("Radius of gaussian filter must be a positive integer.");
3131            }
3132            radius=ex1();
3133            sigma=0.5;
3134            boost::python::extract<double> ex2(filter[2]);
3135            if (!ex2.check() || (sigma=ex2())<=0)
3136            {
3137                throw RipleyException("Sigma must be a positive floating point number.");
3138            }            
3139        }
3140        else
3141        {
3142            throw RipleyException("Unsupported random filter");
3143        }
3144    
3145        // number of points in the internal region
3146        // that is, the ones we need smoothed versions of
3147        const dim_t internal[3] = { m_NN[0], m_NN[1], m_NN[2] };
3148        size_t ext[3];
3149        ext[0]=(size_t)internal[0]+2*radius;  // includes points we need as input
3150        ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3151        ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3152        
3153        // now we check to see if the radius is acceptable
3154        // That is, would not cross multiple ranks in MPI
3155    
3156        if (2*radius>=internal[0]-4)
3157        {
3158            throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
3159        }
3160        if (2*radius>=internal[1]-4)
3161        {
3162            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3163        }
3164        if (2*radius>=internal[2]-4)
3165        {
3166            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3167        }    
3168      
3169        double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3170        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3171        
3172    #ifdef ESYS_MPI
3173      
3174        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3175        {
3176        // since the dimensions are equal for all ranks, this exception
3177        // will be thrown on all ranks
3178        throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3179    
3180        }
3181        dim_t X=m_mpiInfo->rank%m_NX[0];
3182        dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
3183        dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3184    #endif    
3185    
3186    /*    
3187        // if we wanted to test a repeating pattern
3188        size_t basex=0;
3189        size_t basey=0;
3190        size_t basez=0;
3191    #ifdef ESYS_MPI    
3192        basex=X*m_gNE[0]/m_NX[0];
3193        basey=Y*m_gNE[1]/m_NX[1];
3194        basez=Z*m_gNE[2]/m_NX[2];
3195        
3196    cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
3197        
3198    #endif    
3199        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3200    */
3201    
3202    #ifdef ESYS_MPI
3203    
3204    
3205    
3206        BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);
3207        size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3208            // there is an overlap of two points both of which need to have "radius" points on either side.
3209        
3210        size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets
3211        size_t ymidlen=ext[1]-2*inset;  
3212        size_t zmidlen=ext[2]-2*inset;
3213        
3214        Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);    
3215        
3216        MPI_Request reqs[50];       // a non-tight upper bound on how many we need
3217        MPI_Status stats[50];
3218        short rused=0;
3219        
3220        messvec incoms;
3221        messvec outcoms;
3222        
3223        grid.generateInNeighbours(X, Y, Z ,incoms);
3224        grid.generateOutNeighbours(X, Y, Z, outcoms);
3225        
3226        
3227        block.copyAllToBuffer(src);
3228        
3229        int comserr=0;    
3230        for (size_t i=0;i<incoms.size();++i)
3231        {
3232            message& m=incoms[i];
3233            comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3234            block.setUsed(m.destbuffid);
3235        }
3236    
3237        for (size_t i=0;i<outcoms.size();++i)
3238        {
3239            message& m=outcoms[i];
3240            comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3241        }    
3242        
3243        if (!comserr)
3244        {    
3245            comserr=MPI_Waitall(rused, reqs, stats);
3246        }
3247    
3248        if (comserr)
3249        {
3250        // Yes this is throwing an exception as a result of an MPI error.
3251        // and no we don't inform the other ranks that we are doing this.
3252        // however, we have no reason to believe coms work at this point anyway
3253            throw RipleyException("Error in coms for randomFill");      
3254        }
3255        
3256        block.copyUsedFromBuffer(src);
3257        
3258    #endif    
3259        
3260        if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
3261        {
3262          
3263            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3264            escript::Data resdat(0, shape, fs , true);
3265            // don't need to check for exwrite because we just made it
3266            escript::DataVector& dv=resdat.getExpandedVectorReference();
3267        
3268        
3269            // now we need to copy values over
3270            for (size_t z=0;z<(internal[2]);++z)
3271            {
3272                for (size_t y=0;y<(internal[1]);++y)    
3273                {
3274                    for (size_t x=0;x<(internal[0]);++x)
3275                    {
3276                        for (unsigned int i=0;i<numvals;++i)
3277                        {
3278                            dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3279                        }
3280                    }
3281                }
3282            }  
3283            delete[] src;
3284            return resdat;      
3285        }
3286        else        // filter enabled  
3287        {
3288        
3289            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3290            escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3291            // don't need to check for exwrite because we just made it
3292            escript::DataVector& dv=resdat.getExpandedVectorReference();
3293            double* convolution=get3DGauss(radius, sigma);
3294    
3295            for (size_t z=0;z<(internal[2]);++z)
3296            {
3297                for (size_t y=0;y<(internal[1]);++y)    
3298                {
3299                    for (size_t x=0;x<(internal[0]);++x)
3300                    {    
3301                        dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3302                
3303                    }
3304                }
3305            }
3306        
3307            delete[] convolution;
3308            delete[] src;
3309            return resdat;
3310        
3311        }
3312    }
3313    
3314    
3315    
3316    
3317    
3318    
3319  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const {
3320      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3321      //is the found element even owned by this rank      //is the found element even owned by this rank
3322        // (inside owned or shared elements but will map to an owned element)
3323      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3324          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3325                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
3326            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3327                    + m_dx[dim]/2.;
3328            if (min > coords[dim] || max < coords[dim]) {
3329              return NOT_MINE;              return NOT_MINE;
3330          }          }
3331      }      }
# Line 2876  int Brick::findNode(const double *coords Line 3344  int Brick::findNode(const double *coords
3344          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
3345      }      }
3346      //find the closest node      //find the closest node
3347      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
3348          double xdist = x - (ex + dx)*m_dx[0];          double xdist = x - (ex + dx)*m_dx[0];
3349          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
3350              double ydist = y - (ey + dy)*m_dx[1];              double ydist = y - (ey + dy)*m_dx[1];
3351              for (int dz = 0; dz < 2; dz++) {              for (int dz = 0; dz < 1; dz++) {
3352                  double zdist = z - (ez + dz)*m_dx[2];                  double zdist = z - (ez + dz)*m_dx[2];
3353                  double total = xdist*xdist + ydist*ydist + zdist*zdist;                  double total = xdist*xdist + ydist*ydist + zdist*zdist;
3354                  if (total < minDist) {                  if (total < minDist) {
3355                      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],
3356                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3357                      minDist = total;                      minDist = total;
3358                  }                  }
3359              }              }
3360          }          }
3361      }      }
3362        if (closest == NOT_MINE) {
3363            throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");
3364        }
3365      return closest;      return closest;
3366  }  }
3367    
3368  void Brick::setAssembler(std::string type, std::map<std::string,  void Brick::setAssembler(std::string type, std::map<std::string,
3369          escript::Data> constants) {          escript::Data> constants) {
3370      if (type.compare("WaveAssembler") == 0) {      if (type.compare("WaveAssembler") == 0) {
3371            if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3372                throw RipleyException("Domain already using a different custom assembler");
3373            assembler_type = WAVE_ASSEMBLER;
3374          delete assembler;          delete assembler;
3375          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);
3376        } else if (type.compare("LameAssembler") == 0) {
3377            if (assembler_type != LAME_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3378                throw RipleyException("Domain already using a different custom assembler");
3379            assembler_type = LAME_ASSEMBLER;
3380            delete assembler;
3381            assembler = new LameAssembler3D(this, m_dx, m_NX, m_NE, m_NN);
3382      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3383          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Brick does not support the"
3384                                  " requested assembler");                                  " requested assembler");
3385      }      }
3386  }  }

Legend:
Removed from v.4638  
changed lines
  Added in v.4816

  ViewVC Help
Powered by ViewVC 1.1.26