/[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

trunk/ripley/src/Brick.cpp revision 4722 by sshaw, Wed Mar 5 05:29:25 2014 UTC branches/diaplayground/ripley/src/Brick.cpp revision 4988 by caltinay, Wed Jun 4 01:15:10 2014 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
17    #include <limits>
18    
19  #include <ripley/Brick.h>  #include <ripley/Brick.h>
 #include <paso/SystemMatrix.h>  
20  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
21  #include <ripley/DefaultAssembler3D.h>  #include <ripley/DefaultAssembler3D.h>
22  #include <ripley/WaveAssembler3D.h>  #include <ripley/WaveAssembler3D.h>
# Line 60  int indexOfMax(int a, int b, int c) { Line 61  int indexOfMax(int a, int b, int c) {
61  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,
62               double x1, double y1, double z1, int d0, int d1, int d2,               double x1, double y1, double z1, int d0, int d1, int d2,
63               const std::vector<double>& points, const std::vector<int>& tags,               const std::vector<double>& points, const std::vector<int>& tags,
64               const simap_t& tagnamestonums) :               const simap_t& tagnamestonums,
65      RipleyDomain(3)               escript::SubWorld_ptr w) :
66  {      RipleyDomain(3, w)
67    {
68        if (static_cast<long>(n0 + 1) * static_cast<long>(n1 + 1)
69                * static_cast<long>(n2 + 1) > std::numeric_limits<int>::max())
70            throw RipleyException("The number of elements has overflowed, this "
71                    "limit may be raised in future releases.");
72    
73        if (n0 <= 0 || n1 <= 0 || n2 <= 0)
74            throw RipleyException("Number of elements in each spatial dimension "
75                    "must be positive");
76    
77      // ignore subdivision parameters for serial run      // ignore subdivision parameters for serial run
78      if (m_mpiInfo->size == 1) {      if (m_mpiInfo->size == 1) {
79          d0=1;          d0=1;
# Line 192  Brick::Brick(int n0, int n1, int n2, dou Line 203  Brick::Brick(int n0, int n1, int n2, dou
203      populateSampleIds();      populateSampleIds();
204      createPattern();      createPattern();
205            
     assembler = new DefaultAssembler3D(this, m_dx, m_NX, m_NE, m_NN);  
206      for (map<string, int>::const_iterator i = tagnamestonums.begin();      for (map<string, int>::const_iterator i = tagnamestonums.begin();
207              i != tagnamestonums.end(); i++) {              i != tagnamestonums.end(); i++) {
208          setTagMap(i->first, i->second);          setTagMap(i->first, i->second);
# Line 203  Brick::Brick(int n0, int n1, int n2, dou Line 213  Brick::Brick(int n0, int n1, int n2, dou
213    
214  Brick::~Brick()  Brick::~Brick()
215  {  {
     Paso_SystemMatrixPattern_free(m_pattern);  
     Paso_Connector_free(m_connector);  
     delete assembler;  
216  }  }
217    
218  string Brick::getDescription() const  string Brick::getDescription() const
# Line 374  void Brick::readNcGrid(escript::Data& ou Line 381  void Brick::readNcGrid(escript::Data& ou
381  #endif  #endif
382  }  }
383    
384    #ifdef USE_BOOSTIO
385    void Brick::readBinaryGridFromZipped(escript::Data& out, string filename,
386                               const ReaderParameters& params) const
387    {
388        // the mapping is not universally correct but should work on our
389        // supported platforms
390        switch (params.dataType) {
391            case DATATYPE_INT32:
392                readBinaryGridZippedImpl<int>(out, filename, params);
393                break;
394            case DATATYPE_FLOAT32:
395                readBinaryGridZippedImpl<float>(out, filename, params);
396                break;
397            case DATATYPE_FLOAT64:
398                readBinaryGridZippedImpl<double>(out, filename, params);
399                break;
400            default:
401                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
402        }
403    }
404    #endif
405    
406  void Brick::readBinaryGrid(escript::Data& out, string filename,  void Brick::readBinaryGrid(escript::Data& out, string filename,
407                             const ReaderParameters& params) const                             const ReaderParameters& params) const
408  {  {
# Line 423  void Brick::readBinaryGridImpl(escript:: Line 452  void Brick::readBinaryGridImpl(escript::
452      for (size_t i=0; i<params.multiplier.size(); i++)      for (size_t i=0; i<params.multiplier.size(); i++)
453          if (params.multiplier[i]<1)          if (params.multiplier[i]<1)
454              throw RipleyException("readBinaryGrid(): all multipliers must be positive");              throw RipleyException("readBinaryGrid(): all multipliers must be positive");
455        if (params.reverse[0] != 0 || params.reverse[1] != 0)
456            throw RipleyException("readBinaryGrid(): reversing only supported in Z-direction currently");
457    
458      // check file existence and size      // check file existence and size
459      ifstream f(filename.c_str(), ifstream::binary);      ifstream f(filename.c_str(), ifstream::binary);
# Line 455  void Brick::readBinaryGridImpl(escript:: Line 486  void Brick::readBinaryGridImpl(escript::
486      const int first0 = max(0, params.first[0]-m_offset[0]);      const int first0 = max(0, params.first[0]-m_offset[0]);
487      const int first1 = max(0, params.first[1]-m_offset[1]);      const int first1 = max(0, params.first[1]-m_offset[1]);
488      const int first2 = max(0, params.first[2]-m_offset[2]);      const int first2 = max(0, params.first[2]-m_offset[2]);
489        // indices to first value in file (not accounting for reverse yet)
490        int idx0 = max(0, m_offset[0]-params.first[0]);
491        int idx1 = max(0, m_offset[1]-params.first[1]);
492        int idx2 = max(0, m_offset[2]-params.first[2]);
493        // number of values to read
494        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
495        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
496        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
497    
498        // make sure we read the right block if going backwards through file
499        if (params.reverse[2])
500            idx2 = params.numValues[2]-idx2-1;
501    
502        // helpers for reversing
503        const int z_mult = (params.reverse[2] ? -1 : 1);
504    
505        out.requireWrite();
506        vector<ValueType> values(num0*numComp);
507        const int dpp = out.getNumDataPointsPerSample();
508    
509        for (int z=0; z<num2; z++) {
510            for (int y=0; y<num1; y++) {
511                const int fileofs = numComp*(idx0 +
512                                    (idx1+y)*params.numValues[0] +
513                                    (idx2+z_mult*z)*params.numValues[0]*params.numValues[1]);
514                f.seekg(fileofs*sizeof(ValueType));
515                f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
516    
517                for (int x=0; x<num0; x++) {
518                    const int baseIndex = first0+x*params.multiplier[0]
519                                         +(first1+y*params.multiplier[1])*myN0
520                                         +(first2+z*params.multiplier[2])*myN0*myN1;
521                    for (int m2=0; m2<params.multiplier[2]; m2++) {
522                        for (int m1=0; m1<params.multiplier[1]; m1++) {
523                            for (int m0=0; m0<params.multiplier[0]; m0++) {
524                                const int dataIndex = baseIndex+m0
525                                               +m1*myN0
526                                               +m2*myN0*myN1;
527                                double* dest = out.getSampleDataRW(dataIndex);
528                                for (int c=0; c<numComp; c++) {
529                                    ValueType val = values[x*numComp+c];
530    
531                                    if (params.byteOrder != BYTEORDER_NATIVE) {
532                                        char* cval = reinterpret_cast<char*>(&val);
533                                        // this will alter val!!
534                                        byte_swap32(cval);
535                                    }
536                                    if (!std::isnan(val)) {
537                                        for (int q=0; q<dpp; q++) {
538                                            *dest++ = static_cast<double>(val);
539                                        }
540                                    }
541                                }
542                            }
543                        }
544                    }
545                }
546            }
547        }
548    
549        f.close();
550    }
551    
552    #ifdef USE_BOOSTIO
553    template<typename ValueType>
554    void Brick::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
555                                   const ReaderParameters& params) const
556    {
557        // check destination function space
558        int myN0, myN1, myN2;
559        if (out.getFunctionSpace().getTypeCode() == Nodes) {
560            myN0 = m_NN[0];
561            myN1 = m_NN[1];
562            myN2 = m_NN[2];
563        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
564                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
565            myN0 = m_NE[0];
566            myN1 = m_NE[1];
567            myN2 = m_NE[2];
568        } else
569            throw RipleyException("readBinaryGridFromZipped(): invalid function space for output data object");
570    
571        if (params.first.size() != 3)
572            throw RipleyException("readBinaryGridFromZipped(): argument 'first' must have 3 entries");
573    
574        if (params.numValues.size() != 3)
575            throw RipleyException("readBinaryGridFromZipped(): argument 'numValues' must have 3 entries");
576    
577        if (params.multiplier.size() != 3)
578            throw RipleyException("readBinaryGridFromZipped(): argument 'multiplier' must have 3 entries");
579        for (size_t i=0; i<params.multiplier.size(); i++)
580            if (params.multiplier[i]<1)
581                throw RipleyException("readBinaryGridFromZipped(): all multipliers must be positive");
582    
583        // check file existence and size
584        ifstream f(filename.c_str(), ifstream::binary);
585        if (f.fail()) {
586            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
587        }
588        f.seekg(0, ios::end);
589        const int numComp = out.getDataPointSize();
590        int filesize = f.tellg();
591        f.seekg(0, ios::beg);
592        std::vector<char> compressed(filesize);
593        f.read((char*)&compressed[0], filesize);
594        f.close();
595        std::vector<char> decompressed = unzip(compressed);
596        filesize = decompressed.size();
597        const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
598        if (filesize < reqsize) {
599            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
600        }
601    
602        // check if this rank contributes anything
603        if (params.first[0] >= m_offset[0]+myN0 ||
604                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
605                params.first[1] >= m_offset[1]+myN1 ||
606                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
607                params.first[2] >= m_offset[2]+myN2 ||
608                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
609            return;
610        }
611    
612        // now determine how much this rank has to write
613    
614        // first coordinates in data object to write to
615        const int first0 = max(0, params.first[0]-m_offset[0]);
616        const int first1 = max(0, params.first[1]-m_offset[1]);
617        const int first2 = max(0, params.first[2]-m_offset[2]);
618      // indices to first value in file      // indices to first value in file
619      const int idx0 = max(0, m_offset[0]-params.first[0]);      const int idx0 = max(0, m_offset[0]-params.first[0]);
620      const int idx1 = max(0, m_offset[1]-params.first[1]);      const int idx1 = max(0, m_offset[1]-params.first[1]);
# Line 472  void Brick::readBinaryGridImpl(escript:: Line 632  void Brick::readBinaryGridImpl(escript::
632          for (int y=0; y<num1; y++) {          for (int y=0; y<num1; y++) {
633              const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]              const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
634                               +(idx2+z)*params.numValues[0]*params.numValues[1]);                               +(idx2+z)*params.numValues[0]*params.numValues[1]);
635              f.seekg(fileofs*sizeof(ValueType));              memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
636              f.read((char*)&values[0], num0*numComp*sizeof(ValueType));              
   
637              for (int x=0; x<num0; x++) {              for (int x=0; x<num0; x++) {
638                  const int baseIndex = first0+x*params.multiplier[0]                  const int baseIndex = first0+x*params.multiplier[0]
639                                       +(first1+y*params.multiplier[1])*myN0                                       +(first1+y*params.multiplier[1])*myN0
# Line 506  void Brick::readBinaryGridImpl(escript:: Line 665  void Brick::readBinaryGridImpl(escript::
665              }              }
666          }          }
667      }      }
   
     f.close();  
668  }  }
669    #endif
670    
671  void Brick::writeBinaryGrid(const escript::Data& in, string filename,  void Brick::writeBinaryGrid(const escript::Data& in, string filename,
672                              int byteOrder, int dataType) const                              int byteOrder, int dataType) const
# Line 1848  void Brick::assembleIntegrate(vector<dou Line 2006  void Brick::assembleIntegrate(vector<dou
2006  }  }
2007    
2008  //protected  //protected
2009  dim_t Brick::insertNeighbourNodes(IndexVector& index, index_t node) const  IndexVector Brick::getDiagonalIndices() const
2010  {  {
2011        IndexVector ret;
2012      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2013      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
     const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];  
     const int x=node%nDOF0;  
     const int y=node%(nDOF0*nDOF1)/nDOF0;  
     const int z=node/(nDOF0*nDOF1);  
     int num=0;  
     // loop through potential neighbours and add to index if positions are  
     // within bounds  
2014      for (int i2=-1; i2<2; i2++) {      for (int i2=-1; i2<2; i2++) {
2015          for (int i1=-1; i1<2; i1++) {          for (int i1=-1; i1<2; i1++) {
2016              for (int i0=-1; i0<2; i0++) {              for (int i0=-1; i0<2; i0++) {
2017                  // skip node itself                  ret.push_back(i2*nDOF0*nDOF1+i1*nDOF0+i0);
                 if (i0==0 && i1==0 && i2==0)  
                     continue;  
                 // location of neighbour node  
                 const int nx=x+i0;  
                 const int ny=y+i1;  
                 const int nz=z+i2;  
                 if (nx>=0 && ny>=0 && nz>=0  
                         && nx<nDOF0 && ny<nDOF1 && nz<nDOF2) {  
                     index.push_back(nz*nDOF0*nDOF1+ny*nDOF0+nx);  
                     num++;  
                 }  
2018              }              }
2019          }          }
2020      }      }
2021    
2022      return num;      return ret;
2023  }  }
2024    
2025  //protected  //protected
# Line 1909  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_ptr coupler(new paso::Coupler(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));      coupler->startCollect(in.getDataRO());
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 = coupler->finishCollect();
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 1925  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      }      }
     Paso_Coupler_free(coupler);  
2069  }  }
2070    
2071  //private  //private
# Line 1946  void Brick::populateSampleIds() Line 2086  void Brick::populateSampleIds()
2086          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2087      }      }
2088      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2089      m_nodeId.resize(getNumNodes());      
2090      m_dofId.resize(numDOF);      try {
2091      m_elementId.resize(getNumElements());          m_nodeId.resize(getNumNodes());
2092            m_dofId.resize(numDOF);
2093            m_elementId.resize(getNumElements());
2094        } catch (const std::length_error& le) {
2095            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2096        }
2097        
2098      // populate face element counts      // populate face element counts
2099      //left      //left
2100      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2422  void Brick::createPattern() Line 2567  void Brick::createPattern()
2567      }      }
2568    
2569      // create connector      // create connector
2570      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
2571              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2572              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2573      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2574              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2575              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo));
2576      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);  
   
     // create main and couple blocks  
     Paso_Pattern *mainPattern = createMainPattern();  
     Paso_Pattern *colPattern, *rowPattern;  
     createCouplePatterns(colIndices, rowIndices, numShared, &colPattern, &rowPattern);  
   
     // allocate paso distribution  
     Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,  
             const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);  
   
     // finally create the system matrix  
     m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,  
             distribution, distribution, mainPattern, colPattern, rowPattern,  
             m_connector, m_connector);  
   
     Paso_Distribution_free(distribution);  
2577    
2578      // useful debug output      // useful debug output
2579      /*      /*
# Line 2505  void Brick::createPattern() Line 2632  void Brick::createPattern()
2632          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2633      }      }
2634      */      */
   
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colPattern);  
     Paso_Pattern_free(rowPattern);  
2635  }  }
2636    
2637  //private  //private
2638  void Brick::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,  void Brick::addToMatrixAndRHS(SystemMatrix* S, escript::Data& F,
2639           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2640           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2641  {  {
# Line 2536  void Brick::addToMatrixAndRHS(Paso_Syste Line 2659  void Brick::addToMatrixAndRHS(Paso_Syste
2659          }          }
2660      }      }
2661      if (addS) {      if (addS) {
2662          addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);          S->add(rowIndex, EM_S);
2663      }      }
2664  }  }
2665    
# Line 2861  void Brick::interpolateNodesOnFaces(escr Line 2984  void Brick::interpolateNodesOnFaces(escr
2984    
2985  namespace  namespace
2986  {  {
2987      // Calculates a guassian blur colvolution matrix for 3D      // Calculates a gaussian blur convolution matrix for 3D
2988      // See wiki article on the subject      // See wiki article on the subject
2989      double* get3DGauss(unsigned radius, double sigma)      double* get3DGauss(unsigned radius, double sigma)
2990      {      {
# Line 2971  that ripley has. Line 3094  that ripley has.
3094  */  */
3095  escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const  escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
3096  {  {
     if (m_numDim!=3)  
     {  
         throw RipleyException("Brick must be 3D.");  
     }  
       
3097      unsigned int radius=0;  // these are only used by gaussian      unsigned int radius=0;  // these are only used by gaussian
3098      double sigma=0.5;      double sigma=0.5;
3099            
# Line 3002  escript::Data Brick::randomFillWorker(co Line 3120  escript::Data Brick::randomFillWorker(co
3120          boost::python::extract<double> ex2(filter[2]);          boost::python::extract<double> ex2(filter[2]);
3121          if (!ex2.check() || (sigma=ex2())<=0)          if (!ex2.check() || (sigma=ex2())<=0)
3122          {          {
3123              throw RipleyException("Sigma must be a postive floating point number.");              throw RipleyException("Sigma must be a positive floating point number.");
3124          }                      }            
3125      }      }
3126      else      else
3127      {      {
3128          throw RipleyException("Unsupported random filter");          throw RipleyException("Unsupported random filter");
3129      }      }
       
       
3130    
3131            // number of points in the internal region
3132      size_t internal[3];      // that is, the ones we need smoothed versions of
3133      internal[0]=m_NE[0]+1;  // number of points in the internal region      const dim_t internal[3] = { m_NN[0], m_NN[1], m_NN[2] };
     internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of  
     internal[2]=m_NE[2]+1;  // that is, the ones we need smoothed versions of  
3134      size_t ext[3];      size_t ext[3];
3135      ext[0]=internal[0]+2*radius;    // includes points we need as input      ext[0]=(size_t)internal[0]+2*radius;  // includes points we need as input
3136      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3137      ext[2]=internal[2]+2*radius;    // for smoothing      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3138            
3139      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3140      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
# Line 3047  escript::Data Brick::randomFillWorker(co Line 3161  escript::Data Brick::randomFillWorker(co
3161      {      {
3162      // since the dimensions are equal for all ranks, this exception      // since the dimensions are equal for all ranks, this exception
3163      // will be thrown on all ranks      // will be thrown on all ranks
3164      throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");      throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3165    
3166      }      }
3167      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
# Line 3183  cout << "basex=" << basex << " basey=" < Line 3297  cout << "basex=" << basex << " basey=" <
3297      }      }
3298  }  }
3299    
3300    int Brick::findNode(const double *coords) const
3301    {
   
   
   
 int Brick::findNode(const double *coords) const {  
3302      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3303      //is the found element even owned by this rank      //is the found element even owned by this rank
3304      // (inside owned or shared elements but will map to an owned element)      // (inside owned or shared elements but will map to an owned element)
# Line 3205  int Brick::findNode(const double *coords Line 3315  int Brick::findNode(const double *coords
3315      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3316      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3317      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3318        
3319        //check if the point is even inside the domain
3320        if (x < 0 || y < 0 || z < 0
3321                || x > m_length[0] || y > m_length[1] || z > m_length[2])
3322            return NOT_MINE;
3323            
3324      // distance in elements      // distance in elements
3325      int ex = (int) floor(x / m_dx[0]);      int ex = (int) floor(x / m_dx[0]);
3326      int ey = (int) floor(y / m_dx[1]);      int ey = (int) floor(y / m_dx[1]);
# Line 3237  int Brick::findNode(const double *coords Line 3353  int Brick::findNode(const double *coords
3353      return closest;      return closest;
3354  }  }
3355    
3356  void Brick::setAssembler(std::string type, std::map<std::string,  Assembler_ptr Brick::createAssembler(std::string type, std::map<std::string,
3357          escript::Data> constants) {          escript::Data> constants) const
3358      if (type.compare("WaveAssembler") == 0) {  {
3359          delete assembler;      if (type.compare("DefaultAssembler") == 0) {
3360          assembler = new WaveAssembler3D(this, m_dx, m_NX, m_NE, m_NN, constants);          return Assembler_ptr(new DefaultAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
3361        } else if (type.compare("WaveAssembler") == 0) {
3362            return Assembler_ptr(new WaveAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN, constants));
3363        } else if (type.compare("LameAssembler") == 0) {
3364            return Assembler_ptr(new LameAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
3365      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3366          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Brick does not support the"
3367                                  " requested assembler");                                  " requested assembler");
3368      }      }
3369  }  }

Legend:
Removed from v.4722  
changed lines
  Added in v.4988

  ViewVC Help
Powered by ViewVC 1.1.26