/[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 4629 by sshaw, Fri Jan 24 03:29:25 2014 UTC revision 4851 by sshaw, Wed Apr 9 03:30:36 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,
63               const std::map<std::string, int>& 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 51  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 236  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 407  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 543  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 782  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 1940  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      Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));      // expand data object if necessary to be able to grab the whole data
2057        const_cast<escript::Data*>(&in)->expand();
2058        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 1954  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 2222  void Brick::createPattern() Line 2335  void Brick::createPattern()
2335      RankVector neighbour;      RankVector neighbour;
2336      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
2337      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
2338      int numShared=0;      int numShared=0, expectedShared=0;;
2339      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
2340      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];
2341      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 2236  void Brick::createPattern() Line 2349  void Brick::createPattern()
2349                  const int nx=x+i0;                  const int nx=x+i0;
2350                  const int ny=y+i1;                  const int ny=y+i1;
2351                  const int nz=z+i2;                  const int nz=z+i2;
2352                    if (!(nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2])) {
2353                        continue;
2354                    }
2355                    if (i0==0 && i1==0)
2356                        expectedShared += nDOF0*nDOF1;
2357                    else if (i0==0 && i2==0)
2358                        expectedShared += nDOF0*nDOF2;
2359                    else if (i1==0 && i2==0)
2360                        expectedShared += nDOF1*nDOF2;
2361                    else if (i0==0)
2362                        expectedShared += nDOF0;
2363                    else if (i1==0)
2364                        expectedShared += nDOF1;
2365                    else if (i2==0)
2366                        expectedShared += nDOF2;
2367                    else
2368                        expectedShared++;
2369                }
2370            }
2371        }
2372        
2373        vector<IndexVector> rowIndices(expectedShared);
2374        
2375        for (int i2=-1; i2<2; i2++) {
2376            for (int i1=-1; i1<2; i1++) {
2377                for (int i0=-1; i0<2; i0++) {
2378                    // skip this rank
2379                    if (i0==0 && i1==0 && i2==0)
2380                        continue;
2381                    // location of neighbour rank
2382                    const int nx=x+i0;
2383                    const int ny=y+i1;
2384                    const int nz=z+i2;
2385                  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]) {
2386                      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);
2387                      if (i0==0 && i1==0) {                      if (i0==0 && i1==0) {
# Line 2251  void Brick::createPattern() Line 2397  void Brick::createPattern()
2397                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2398                                  if (j>0) {                                  if (j>0) {
2399                                      if (i>0)                                      if (i>0)
2400                                          colIndices[firstDOF+j-1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0, numShared);
2401                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2402                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2403                                          colIndices[firstDOF+j-1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0, numShared);
2404                                  }                                  }
2405                                  if (i>0)                                  if (i>0)
2406                                      colIndices[firstDOF+j-nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0, numShared);
2407                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2408                                  if (i<nDOF1-1)                                  if (i<nDOF1-1)
2409                                      colIndices[firstDOF+j+nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0, numShared);
2410                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2411                                      if (i>0)                                      if (i>0)
2412                                          colIndices[firstDOF+j+1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0, numShared);
2413                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2414                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2415                                          colIndices[firstDOF+j+1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0, numShared);
2416                                  }                                  }
2417                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2418                              }                              }
# Line 2285  void Brick::createPattern() Line 2431  void Brick::createPattern()
2431                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2432                                  if (j>0) {                                  if (j>0) {
2433                                      if (i>0)                                      if (i>0)
2434                                          colIndices[firstDOF+j-1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0*nDOF1, numShared);
2435                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2436                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2437                                          colIndices[firstDOF+j-1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0*nDOF1, numShared);
2438                                  }                                  }
2439                                  if (i>0)                                  if (i>0)
2440                                      colIndices[firstDOF+j-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0*nDOF1, numShared);
2441                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2442                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2443                                      colIndices[firstDOF+j+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0*nDOF1, numShared);
2444                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2445                                      if (i>0)                                      if (i>0)
2446                                          colIndices[firstDOF+j+1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0*nDOF1, numShared);
2447                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2448                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2449                                          colIndices[firstDOF+j+1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0*nDOF1, numShared);
2450                                  }                                  }
2451                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2452                              }                              }
# Line 2319  void Brick::createPattern() Line 2465  void Brick::createPattern()
2465                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2466                                  if (j>0) {                                  if (j>0) {
2467                                      if (i>0)                                      if (i>0)
2468                                          colIndices[firstDOF+(j-1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0-nDOF0*nDOF1, numShared);
2469                                      colIndices[firstDOF+(j-1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0, numShared);
2470                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2471                                          colIndices[firstDOF+(j-1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0+nDOF0*nDOF1, numShared);
2472                                  }                                  }
2473                                  if (i>0)                                  if (i>0)
2474                                      colIndices[firstDOF+j*nDOF0-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0-nDOF0*nDOF1, numShared);
2475                                  colIndices[firstDOF+j*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0, numShared);
2476                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2477                                      colIndices[firstDOF+j*nDOF0+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0+nDOF0*nDOF1, numShared);
2478                                  if (j<nDOF1-1) {                                  if (j<nDOF1-1) {
2479                                      if (i>0)                                      if (i>0)
2480                                          colIndices[firstDOF+(j+1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0-nDOF0*nDOF1, numShared);
2481                                      colIndices[firstDOF+(j+1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0, numShared);
2482                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2483                                          colIndices[firstDOF+(j+1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0+nDOF0*nDOF1, numShared);
2484                                  }                                  }
2485                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2486                              }                              }
# Line 2350  void Brick::createPattern() Line 2496  void Brick::createPattern()
2496                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
2497                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2498                              if (i>0)                              if (i>0)
2499                                  colIndices[firstDOF+i-1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
2500                              colIndices[firstDOF+i].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
2501                              if (i<nDOF0-1)                              if (i<nDOF0-1)
2502                                  colIndices[firstDOF+i+1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
2503                              m_dofMap[firstNode+i]=numDOF+numShared;                              m_dofMap[firstNode+i]=numDOF+numShared;
2504                          }                          }
2505                      } else if (i1==0) {                      } else if (i1==0) {
# Line 2368  void Brick::createPattern() Line 2514  void Brick::createPattern()
2514                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
2515                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2516                              if (i>0)                              if (i>0)
2517                                  colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
2518                              colIndices[firstDOF+i*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
2519                              if (i<nDOF1-1)                              if (i<nDOF1-1)
2520                                  colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
2521                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2522                          }                          }
2523                      } else if (i2==0) {                      } else if (i2==0) {
# Line 2386  void Brick::createPattern() Line 2532  void Brick::createPattern()
2532                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2533                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2534                              if (i>0)                              if (i>0)
2535                                  colIndices[firstDOF+(i-1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0*nDOF1, numShared);
2536                              colIndices[firstDOF+i*nDOF0*nDOF1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0*nDOF1, numShared);
2537                              if (i<nDOF2-1)                              if (i<nDOF2-1)
2538                                  colIndices[firstDOF+(i+1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0*nDOF1, numShared);
2539                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2540                          }                          }
2541                      } else {                      } else {
# Line 2403  void Brick::createPattern() Line 2549  void Brick::createPattern()
2549                          offsetInShared.push_back(offsetInShared.back()+1);                          offsetInShared.push_back(offsetInShared.back()+1);
2550                          sendShared.push_back(dof);                          sendShared.push_back(dof);
2551                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
2552                          colIndices[dof].push_back(numShared);                          doublyLink(colIndices, rowIndices, dof, numShared);
2553                          m_dofMap[node]=numDOF+numShared;                          m_dofMap[node]=numDOF+numShared;
2554                          ++numShared;                          ++numShared;
2555                      }                      }
# Line 2412  void Brick::createPattern() Line 2558  void Brick::createPattern()
2558          }          }
2559      }      }
2560    
2561    #pragma omp parallel for
2562        for (int i = 0; i < numShared; i++) {
2563            std::sort(rowIndices[i].begin(), rowIndices[i].end());
2564        }
2565    
2566      // create connector      // create connector
2567      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
2568              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2569              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2570      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2571              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2572              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2573      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);  
2574    
2575      // create main and couple blocks      // create main and couple blocks
2576      Paso_Pattern *mainPattern = createMainPattern();      paso::Pattern_ptr mainPattern = createMainPattern();
2577      Paso_Pattern *colPattern, *rowPattern;      paso::Pattern_ptr colPattern, rowPattern;
2578      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, colPattern, rowPattern);
2579    
2580      // allocate paso distribution      // allocate paso distribution
2581      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2582              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2583    
2584      // finally create the system matrix      // finally create the system matrix
2585      m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,      m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2586              distribution, distribution, mainPattern, colPattern, rowPattern,              distribution, distribution, mainPattern, colPattern, rowPattern,
2587              m_connector, m_connector);              m_connector, m_connector));
   
     Paso_Distribution_free(distribution);  
2588    
2589      // useful debug output      // useful debug output
2590      /*      /*
# Line 2496  void Brick::createPattern() Line 2643  void Brick::createPattern()
2643          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2644      }      }
2645      */      */
   
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colPattern);  
     Paso_Pattern_free(rowPattern);  
2646  }  }
2647    
2648  //private  //private
2649  void Brick::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,  void Brick::addToMatrixAndRHS(paso::SystemMatrix_ptr S, escript::Data& F,
2650           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2651           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2652  {  {
# Line 2850  void Brick::interpolateNodesOnFaces(escr Line 2993  void Brick::interpolateNodesOnFaces(escr
2993      }      }
2994  }  }
2995    
2996    namespace
2997    {
2998        // Calculates a gaussian blur convolution matrix for 3D
2999        // See wiki article on the subject
3000        double* get3DGauss(unsigned radius, double sigma)
3001        {
3002            double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];
3003            double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);
3004            double total=0;
3005            int r=static_cast<int>(radius);
3006            for (int z=-r;z<=r;++z)
3007            {
3008                for (int y=-r;y<=r;++y)
3009                {
3010                    for (int x=-r;x<=r;++x)
3011                    {        
3012                        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));
3013                        total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];    
3014                    }
3015                }
3016            }
3017            double invtotal=1/total;
3018            for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
3019            {
3020                arr[p]*=invtotal;
3021            }
3022            return arr;
3023        }
3024        
3025        // applies conv to source to get a point.
3026        // (xp, yp) are the coords in the source matrix not the destination matrix
3027        double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)
3028        {
3029            size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3030            size_t sbase=bx+by*width+bz*width*height;
3031            double result=0;
3032            for (int z=0;z<2*radius+1;++z)
3033            {
3034                for (int y=0;y<2*radius+1;++y)
3035                {    
3036                    for (int x=0;x<2*radius+1;++x)
3037                    {
3038                        result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];
3039                    }
3040                }
3041            }
3042            // use this line for "pass-though" (return the centre point value)
3043    //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3044            return result;      
3045        }
3046    }
3047    
3048    /* This is a wrapper for filtered (and non-filtered) randoms
3049     * For detailed doco see randomFillWorker
3050    */
3051    escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
3052           const escript::FunctionSpace& what,
3053           long seed, const boost::python::tuple& filter) const
3054    {
3055        int numvals=escript::DataTypes::noValues(shape);
3056        if (len(filter)>0 && (numvals!=1))
3057        {
3058            throw RipleyException("Ripley only supports filters for scalar data.");
3059        }
3060        escript::Data res=randomFillWorker(shape, seed, filter);
3061        if (res.getFunctionSpace()!=what)
3062        {
3063            escript::Data r=escript::Data(res, what);
3064            return r;
3065        }
3066        return res;
3067    }
3068    
3069    /* This routine produces a Data object filled with smoothed random data.
3070    The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.
3071    A parameter radius  gives the size of the stencil used for the smoothing.
3072    A point on the left hand edge for example, will still require `radius` extra points to the left
3073    in order to complete the stencil.
3074    
3075    All local calculation is done on an array called `src`, with
3076    dimensions = ext[0] * ext[1] *ext[2].
3077    Where ext[i]= internal[i]+2*radius.
3078    
3079    Now for MPI there is overlap to deal with. We need to share both the overlapping
3080    values themselves but also the external region.
3081    
3082    In a hypothetical 1-D case:
3083    
3084    
3085    1234567
3086    would be split into two ranks thus:
3087    123(4)  (4)567     [4 being a shared element]
3088    
3089    If the radius is 2.   There will be padding elements on the outside:
3090    
3091    pp123(4)  (4)567pp
3092    
3093    To ensure that 4 can be correctly computed on both ranks, values from the other rank
3094    need to be known.
3095    
3096    pp123(4)56   23(4)567pp
3097    
3098    Now in our case, we wout set all the values 23456 on the left rank and send them to the
3099    right hand rank.
3100    
3101    So the edges _may_ need to be shared at a distance `inset` from all boundaries.
3102    inset=2*radius+1    
3103    This is to ensure that values at distance `radius` from the shared/overlapped element
3104    that ripley has.
3105    */
3106    escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
3107    {
3108        unsigned int radius=0;  // these are only used by gaussian
3109        double sigma=0.5;
3110        
3111        unsigned int numvals=escript::DataTypes::noValues(shape);
3112        
3113        if (len(filter)==0)
3114        {
3115        // nothing special required here yet
3116        }
3117        else if (len(filter)==3)
3118        {
3119            boost::python::extract<string> ex(filter[0]);
3120            if (!ex.check() || (ex()!="gaussian"))
3121            {
3122                throw RipleyException("Unsupported random filter for Brick.");
3123            }
3124            boost::python::extract<unsigned int> ex1(filter[1]);
3125            if (!ex1.check())
3126            {
3127                throw RipleyException("Radius of gaussian filter must be a positive integer.");
3128            }
3129            radius=ex1();
3130            sigma=0.5;
3131            boost::python::extract<double> ex2(filter[2]);
3132            if (!ex2.check() || (sigma=ex2())<=0)
3133            {
3134                throw RipleyException("Sigma must be a positive floating point number.");
3135            }            
3136        }
3137        else
3138        {
3139            throw RipleyException("Unsupported random filter");
3140        }
3141    
3142        // number of points in the internal region
3143        // that is, the ones we need smoothed versions of
3144        const dim_t internal[3] = { m_NN[0], m_NN[1], m_NN[2] };
3145        size_t ext[3];
3146        ext[0]=(size_t)internal[0]+2*radius;  // includes points we need as input
3147        ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3148        ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3149        
3150        // now we check to see if the radius is acceptable
3151        // That is, would not cross multiple ranks in MPI
3152    
3153        if (2*radius>=internal[0]-4)
3154        {
3155            throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
3156        }
3157        if (2*radius>=internal[1]-4)
3158        {
3159            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3160        }
3161        if (2*radius>=internal[2]-4)
3162        {
3163            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3164        }    
3165      
3166        double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3167        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3168        
3169    #ifdef ESYS_MPI
3170      
3171        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3172        {
3173        // since the dimensions are equal for all ranks, this exception
3174        // will be thrown on all ranks
3175        throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3176    
3177        }
3178        dim_t X=m_mpiInfo->rank%m_NX[0];
3179        dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
3180        dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3181    #endif    
3182    
3183    /*    
3184        // if we wanted to test a repeating pattern
3185        size_t basex=0;
3186        size_t basey=0;
3187        size_t basez=0;
3188    #ifdef ESYS_MPI    
3189        basex=X*m_gNE[0]/m_NX[0];
3190        basey=Y*m_gNE[1]/m_NX[1];
3191        basez=Z*m_gNE[2]/m_NX[2];
3192        
3193    cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
3194        
3195    #endif    
3196        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3197    */
3198    
3199    #ifdef ESYS_MPI
3200    
3201    
3202    
3203        BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);
3204        size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3205            // there is an overlap of two points both of which need to have "radius" points on either side.
3206        
3207        size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets
3208        size_t ymidlen=ext[1]-2*inset;  
3209        size_t zmidlen=ext[2]-2*inset;
3210        
3211        Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);    
3212        
3213        MPI_Request reqs[50];       // a non-tight upper bound on how many we need
3214        MPI_Status stats[50];
3215        short rused=0;
3216        
3217        messvec incoms;
3218        messvec outcoms;
3219        
3220        grid.generateInNeighbours(X, Y, Z ,incoms);
3221        grid.generateOutNeighbours(X, Y, Z, outcoms);
3222        
3223        
3224        block.copyAllToBuffer(src);
3225        
3226        int comserr=0;    
3227        for (size_t i=0;i<incoms.size();++i)
3228        {
3229            message& m=incoms[i];
3230            comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3231            block.setUsed(m.destbuffid);
3232        }
3233    
3234        for (size_t i=0;i<outcoms.size();++i)
3235        {
3236            message& m=outcoms[i];
3237            comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3238        }    
3239        
3240        if (!comserr)
3241        {    
3242            comserr=MPI_Waitall(rused, reqs, stats);
3243        }
3244    
3245        if (comserr)
3246        {
3247        // Yes this is throwing an exception as a result of an MPI error.
3248        // and no we don't inform the other ranks that we are doing this.
3249        // however, we have no reason to believe coms work at this point anyway
3250            throw RipleyException("Error in coms for randomFill");      
3251        }
3252        
3253        block.copyUsedFromBuffer(src);
3254        
3255    #endif    
3256        
3257        if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
3258        {
3259          
3260            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3261            escript::Data resdat(0, shape, fs , true);
3262            // don't need to check for exwrite because we just made it
3263            escript::DataVector& dv=resdat.getExpandedVectorReference();
3264        
3265        
3266            // now we need to copy values over
3267            for (size_t z=0;z<(internal[2]);++z)
3268            {
3269                for (size_t y=0;y<(internal[1]);++y)    
3270                {
3271                    for (size_t x=0;x<(internal[0]);++x)
3272                    {
3273                        for (unsigned int i=0;i<numvals;++i)
3274                        {
3275                            dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3276                        }
3277                    }
3278                }
3279            }  
3280            delete[] src;
3281            return resdat;      
3282        }
3283        else        // filter enabled  
3284        {
3285        
3286            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3287            escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3288            // don't need to check for exwrite because we just made it
3289            escript::DataVector& dv=resdat.getExpandedVectorReference();
3290            double* convolution=get3DGauss(radius, sigma);
3291    
3292            for (size_t z=0;z<(internal[2]);++z)
3293            {
3294                for (size_t y=0;y<(internal[1]);++y)    
3295                {
3296                    for (size_t x=0;x<(internal[0]);++x)
3297                    {    
3298                        dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3299                
3300                    }
3301                }
3302            }
3303        
3304            delete[] convolution;
3305            delete[] src;
3306            return resdat;
3307        
3308        }
3309    }
3310    
3311    
3312    
3313    
3314    
3315    
3316  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const {
3317      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3318      //is the found element even owned by this rank      //is the found element even owned by this rank
3319        // (inside owned or shared elements but will map to an owned element)
3320      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3321          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3322                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
3323            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3324                    + m_dx[dim]/2.;
3325            if (min > coords[dim] || max < coords[dim]) {
3326              return NOT_MINE;              return NOT_MINE;
3327          }          }
3328      }      }
# Line 2863  int Brick::findNode(const double *coords Line 3330  int Brick::findNode(const double *coords
3330      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3331      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3332      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3333        
3334        //check if the point is even inside the domain
3335        if (x < 0 || y < 0 || z < 0
3336                || x > m_length[0] || y > m_length[1] || z > m_length[2])
3337            return NOT_MINE;
3338            
3339      // distance in elements      // distance in elements
3340      int ex = (int) floor(x / m_dx[0]);      int ex = (int) floor(x / m_dx[0]);
3341      int ey = (int) floor(y / m_dx[1]);      int ey = (int) floor(y / m_dx[1]);
# Line 2874  int Brick::findNode(const double *coords Line 3347  int Brick::findNode(const double *coords
3347          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
3348      }      }
3349      //find the closest node      //find the closest node
3350      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
3351          double xdist = x - (ex + dx)*m_dx[0];          double xdist = x - (ex + dx)*m_dx[0];
3352          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
3353              double ydist = y - (ey + dy)*m_dx[1];              double ydist = y - (ey + dy)*m_dx[1];
3354              for (int dz = 0; dz < 2; dz++) {              for (int dz = 0; dz < 1; dz++) {
3355                  double zdist = z - (ez + dz)*m_dx[2];                  double zdist = z - (ez + dz)*m_dx[2];
3356                  double total = xdist*xdist + ydist*ydist + zdist*zdist;                  double total = xdist*xdist + ydist*ydist + zdist*zdist;
3357                  if (total < minDist) {                  if (total < minDist) {
3358                      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],
3359                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3360                      minDist = total;                      minDist = total;
3361                  }                  }
3362              }              }
3363          }          }
3364      }      }
3365        if (closest == NOT_MINE) {
3366            throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");
3367        }
3368      return closest;      return closest;
3369  }  }
3370    
3371  void Brick::setAssembler(std::string type, std::map<std::string,  void Brick::setAssembler(std::string type, std::map<std::string,
3372          escript::Data> constants) {          escript::Data> constants) {
3373      if (type.compare("WaveAssembler") == 0) {      if (type.compare("WaveAssembler") == 0) {
3374            if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3375                throw RipleyException("Domain already using a different custom assembler");
3376            assembler_type = WAVE_ASSEMBLER;
3377          delete assembler;          delete assembler;
3378          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);
3379        } else if (type.compare("LameAssembler") == 0) {
3380            if (assembler_type != LAME_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3381                throw RipleyException("Domain already using a different custom assembler");
3382            assembler_type = LAME_ASSEMBLER;
3383            delete assembler;
3384            assembler = new LameAssembler3D(this, m_dx, m_NX, m_NE, m_NN);
3385      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3386          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Brick does not support the"
3387                                  " requested assembler");                                  " requested assembler");
3388      }      }
3389  }  }

Legend:
Removed from v.4629  
changed lines
  Added in v.4851

  ViewVC Help
Powered by ViewVC 1.1.26