/[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 4819 by caltinay, Tue Apr 1 03:50:23 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  {  {
     delete assembler;  
216  }  }
217    
218  string Brick::getDescription() const  string Brick::getDescription() const
# Line 443  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 475  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      // indices to first value in file (not accounting for reverse yet)
490      const int idx0 = max(0, m_offset[0]-params.first[0]);      int idx0 = max(0, m_offset[0]-params.first[0]);
491      const int idx1 = max(0, m_offset[1]-params.first[1]);      int idx1 = max(0, m_offset[1]-params.first[1]);
492      const int idx2 = max(0, m_offset[2]-params.first[2]);      int idx2 = max(0, m_offset[2]-params.first[2]);
493      // number of values to read      // number of values to read
494      const int num0 = min(params.numValues[0]-idx0, myN0-first0);      const int num0 = min(params.numValues[0]-idx0, myN0-first0);
495      const int num1 = min(params.numValues[1]-idx1, myN1-first1);      const int num1 = min(params.numValues[1]-idx1, myN1-first1);
496      const int num2 = min(params.numValues[2]-idx2, myN2-first2);      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();      out.requireWrite();
506      vector<ValueType> values(num0*numComp);      vector<ValueType> values(num0*numComp);
507      const int dpp = out.getNumDataPointsPerSample();      const int dpp = out.getNumDataPointsPerSample();
508    
509      for (int z=0; z<num2; z++) {      for (int z=0; z<num2; z++) {
510          for (int y=0; y<num1; y++) {          for (int y=0; y<num1; y++) {
511              const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]              const int fileofs = numComp*(idx0 +
512                               +(idx2+z)*params.numValues[0]*params.numValues[1]);                                  (idx1+y)*params.numValues[0] +
513                                    (idx2+z_mult*z)*params.numValues[0]*params.numValues[1]);
514              f.seekg(fileofs*sizeof(ValueType));              f.seekg(fileofs*sizeof(ValueType));
515              f.read((char*)&values[0], num0*numComp*sizeof(ValueType));              f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
516    
# Line 1987  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 2051  void Brick::dofToNodes(escript::Data& ou Line 2053  void Brick::dofToNodes(escript::Data& ou
2053      paso::Coupler_ptr coupler(new paso::Coupler(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      coupler->startCollect(in.getSampleDataRO(0));      coupler->startCollect(in.getDataRO());
2057    
2058      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2059      out.requireWrite();      out.requireWrite();
# Line 2084  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 2568  void Brick::createPattern() Line 2575  void Brick::createPattern()
2575              &offsetInShared[0], 1, 0, m_mpiInfo));              &offsetInShared[0], 1, 0, m_mpiInfo));
2576      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
2577    
     // create main and couple blocks  
     paso::Pattern_ptr mainPattern = createMainPattern();  
     paso::Pattern_ptr colPattern, rowPattern;  
     createCouplePatterns(colIndices, rowIndices, numShared, colPattern, rowPattern);  
   
     // allocate paso distribution  
     paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,  
             const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));  
   
     // finally create the system matrix  
     m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,  
             distribution, distribution, mainPattern, colPattern, rowPattern,  
             m_connector, m_connector));  
   
2578      // useful debug output      // useful debug output
2579      /*      /*
2580      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
# Line 2642  void Brick::createPattern() Line 2635  void Brick::createPattern()
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 2666  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 3304  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 3326  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 3358  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          if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)      if (type.compare("DefaultAssembler") == 0) {
3360              throw RipleyException("Domain already using a different custom assembler");          return Assembler_ptr(new DefaultAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
3361          assembler_type = WAVE_ASSEMBLER;      } else if (type.compare("WaveAssembler") == 0) {
3362          delete assembler;          return Assembler_ptr(new WaveAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN, constants));
         assembler = new WaveAssembler3D(this, m_dx, m_NX, m_NE, m_NN, constants);  
3363      } else if (type.compare("LameAssembler") == 0) {      } else if (type.compare("LameAssembler") == 0) {
3364          if (assembler_type != LAME_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)          return Assembler_ptr(new LameAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
             throw RipleyException("Domain already using a different custom assembler");  
         assembler_type = LAME_ASSEMBLER;  
         delete assembler;  
         assembler = new LameAssembler3D(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::Brick does not support the"          throw RipleyException("Ripley::Brick does not support the"
3367                                  " requested assembler");                                  " requested assembler");

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

  ViewVC Help
Powered by ViewVC 1.1.26