/[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 4940 by caltinay, Thu May 15 01:40:06 2014 UTC revision 4988 by caltinay, Wed Jun 4 01:15:10 2014 UTC
# Line 17  Line 17 
17  #include <limits>  #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 204  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 215  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 455  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 487  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 1999  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 2585  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 2659  void Brick::createPattern() Line 2635  void Brick::createPattern()
2635  }  }
2636    
2637  //private  //private
2638  void Brick::addToMatrixAndRHS(paso::SystemMatrix_ptr 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 2683  void Brick::addToMatrixAndRHS(paso::Syst Line 2659  void Brick::addToMatrixAndRHS(paso::Syst
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 3321  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 3381  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.4940  
changed lines
  Added in v.4988

  ViewVC Help
Powered by ViewVC 1.1.26