/[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 4941 by caltinay, Thu May 15 01:49:48 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 1991  void Brick::assembleIntegrate(vector<dou Line 1998  void Brick::assembleIntegrate(vector<dou
1998  }  }
1999    
2000  //protected  //protected
2001  dim_t Brick::insertNeighbourNodes(IndexVector& index, index_t node) const  IndexVector Brick::getDiagonalIndices() const
2002  {  {
2003        IndexVector ret;
2004      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2005      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  
2006      for (int i2=-1; i2<2; i2++) {      for (int i2=-1; i2<2; i2++) {
2007          for (int i1=-1; i1<2; i1++) {          for (int i1=-1; i1<2; i1++) {
2008              for (int i0=-1; i0<2; i0++) {              for (int i0=-1; i0<2; i0++) {
2009                  // 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++;  
                 }  
2010              }              }
2011          }          }
2012      }      }
2013    
2014      return num;      return ret;
2015  }  }
2016    
2017  //protected  //protected
# Line 2055  void Brick::dofToNodes(escript::Data& ou Line 2045  void Brick::dofToNodes(escript::Data& ou
2045      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2046      // 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
2047      const_cast<escript::Data*>(&in)->expand();      const_cast<escript::Data*>(&in)->expand();
2048      coupler->startCollect(in.getSampleDataRO(0));      coupler->startCollect(in.getDataRO());
2049    
2050      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2051      out.requireWrite();      out.requireWrite();
# Line 2088  void Brick::populateSampleIds() Line 2078  void Brick::populateSampleIds()
2078          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2079      }      }
2080      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2081      m_nodeId.resize(getNumNodes());      
2082      m_dofId.resize(numDOF);      try {
2083      m_elementId.resize(getNumElements());          m_nodeId.resize(getNumNodes());
2084            m_dofId.resize(numDOF);
2085            m_elementId.resize(getNumElements());
2086        } catch (const std::length_error& le) {
2087            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2088        }
2089        
2090      // populate face element counts      // populate face element counts
2091      //left      //left
2092      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2572  void Brick::createPattern() Line 2567  void Brick::createPattern()
2567              &offsetInShared[0], 1, 0, m_mpiInfo));              &offsetInShared[0], 1, 0, m_mpiInfo));
2568      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
2569    
     // 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));  
   
2570      // useful debug output      // useful debug output
2571      /*      /*
2572      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
# Line 2646  void Brick::createPattern() Line 2627  void Brick::createPattern()
2627  }  }
2628    
2629  //private  //private
2630  void Brick::addToMatrixAndRHS(paso::SystemMatrix_ptr S, escript::Data& F,  void Brick::addToMatrixAndRHS(SystemMatrix* S, escript::Data& F,
2631           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2632           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2633  {  {
# Line 2670  void Brick::addToMatrixAndRHS(paso::Syst Line 2651  void Brick::addToMatrixAndRHS(paso::Syst
2651          }          }
2652      }      }
2653      if (addS) {      if (addS) {
2654          addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);          S->add(rowIndex, EM_S);
2655      }      }
2656  }  }
2657    

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

  ViewVC Help
Powered by ViewVC 1.1.26