/[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 4851 by sshaw, Wed Apr 9 03:30:36 2014 UTC branches/diaplayground/ripley/src/Brick.cpp revision 4949 by caltinay, Mon May 19 05:54:58 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)      if (n0 <= 0 || n1 <= 0 || n2 <= 0)
74          throw RipleyException("Number of elements in each spatial dimension "          throw RipleyException("Number of elements in each spatial dimension "
75                  "must be positive");                  "must be positive");
# Line 196  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 207  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 1991  void Brick::assembleIntegrate(vector<dou Line 1996  void Brick::assembleIntegrate(vector<dou
1996  }  }
1997    
1998  //protected  //protected
1999  dim_t Brick::insertNeighbourNodes(IndexVector& index, index_t node) const  IndexVector Brick::getDiagonalIndices() const
2000  {  {
2001        IndexVector ret;
2002      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2003      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  
2004      for (int i2=-1; i2<2; i2++) {      for (int i2=-1; i2<2; i2++) {
2005          for (int i1=-1; i1<2; i1++) {          for (int i1=-1; i1<2; i1++) {
2006              for (int i0=-1; i0<2; i0++) {              for (int i0=-1; i0<2; i0++) {
2007                  // 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++;  
                 }  
2008              }              }
2009          }          }
2010      }      }
2011    
2012      return num;      return ret;
2013  }  }
2014    
2015  //protected  //protected
# Line 2055  void Brick::dofToNodes(escript::Data& ou Line 2043  void Brick::dofToNodes(escript::Data& ou
2043      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2044      // 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
2045      const_cast<escript::Data*>(&in)->expand();      const_cast<escript::Data*>(&in)->expand();
2046      coupler->startCollect(in.getSampleDataRO(0));      coupler->startCollect(in.getDataRO());
2047    
2048      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2049      out.requireWrite();      out.requireWrite();
# Line 2088  void Brick::populateSampleIds() Line 2076  void Brick::populateSampleIds()
2076          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2077      }      }
2078      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2079      m_nodeId.resize(getNumNodes());      
2080      m_dofId.resize(numDOF);      try {
2081      m_elementId.resize(getNumElements());          m_nodeId.resize(getNumNodes());
2082            m_dofId.resize(numDOF);
2083            m_elementId.resize(getNumElements());
2084        } catch (const std::length_error& le) {
2085            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2086        }
2087        
2088      // populate face element counts      // populate face element counts
2089      //left      //left
2090      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2572  void Brick::createPattern() Line 2565  void Brick::createPattern()
2565              &offsetInShared[0], 1, 0, m_mpiInfo));              &offsetInShared[0], 1, 0, m_mpiInfo));
2566      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
2567    
     // 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));  
   
2568      // useful debug output      // useful debug output
2569      /*      /*
2570      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
# Line 2646  void Brick::createPattern() Line 2625  void Brick::createPattern()
2625  }  }
2626    
2627  //private  //private
2628  void Brick::addToMatrixAndRHS(paso::SystemMatrix_ptr S, escript::Data& F,  void Brick::addToMatrixAndRHS(SystemMatrix* S, escript::Data& F,
2629           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2630           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2631  {  {
# Line 2670  void Brick::addToMatrixAndRHS(paso::Syst Line 2649  void Brick::addToMatrixAndRHS(paso::Syst
2649          }          }
2650      }      }
2651      if (addS) {      if (addS) {
2652          addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);          S->add(rowIndex, EM_S);
2653      }      }
2654  }  }
2655    
# Line 3308  cout << "basex=" << basex << " basey=" < Line 3287  cout << "basex=" << basex << " basey=" <
3287      }      }
3288  }  }
3289    
3290    int Brick::findNode(const double *coords) const
3291    {
   
   
   
 int Brick::findNode(const double *coords) const {  
3292      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3293      //is the found element even owned by this rank      //is the found element even owned by this rank
3294      // (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 3368  int Brick::findNode(const double *coords Line 3343  int Brick::findNode(const double *coords
3343      return closest;      return closest;
3344  }  }
3345    
3346  void Brick::setAssembler(std::string type, std::map<std::string,  Assembler_ptr Brick::createAssembler(std::string type, std::map<std::string,
3347          escript::Data> constants) {          escript::Data> constants) const
3348      if (type.compare("WaveAssembler") == 0) {  {
3349          if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)      if (type.compare("DefaultAssembler") == 0) {
3350              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));
3351          assembler_type = WAVE_ASSEMBLER;      } else if (type.compare("WaveAssembler") == 0) {
3352          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);  
3353      } else if (type.compare("LameAssembler") == 0) {      } else if (type.compare("LameAssembler") == 0) {
3354          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);  
3355      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3356          throw RipleyException("Ripley::Brick does not support the"          throw RipleyException("Ripley::Brick does not support the"
3357                                  " requested assembler");                                  " requested assembler");

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

  ViewVC Help
Powered by ViewVC 1.1.26