/[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 4816 by caltinay, Fri Mar 28 06:16:02 2014 UTC revision 4861 by sshaw, Thu Apr 10 05:17:47 2014 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
17    #include <limits>
18    
19  #include <ripley/Brick.h>  #include <ripley/Brick.h>
20  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
21  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
# Line 63  Brick::Brick(int n0, int n1, int n2, dou Line 65  Brick::Brick(int n0, int n1, int n2, dou
65               const simap_t& tagnamestonums) :               const simap_t& tagnamestonums) :
66      RipleyDomain(3)      RipleyDomain(3)
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 203  Brick::Brick(int n0, int n1, int n2, dou Line 214  Brick::Brick(int n0, int n1, int n2, dou
214    
215  Brick::~Brick()  Brick::~Brick()
216  {  {
     paso::SystemMatrixPattern_free(m_pattern);  
     paso::Connector_free(m_connector);  
217      delete assembler;      delete assembler;
218  }  }
219    
# Line 2050  void Brick::nodesToDOF(escript::Data& ou Line 2059  void Brick::nodesToDOF(escript::Data& ou
2059  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const
2060  {  {
2061      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
2062      paso::Coupler* coupler = paso::Coupler_alloc(m_connector, numComp);      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2063      // 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
2064      const_cast<escript::Data*>(&in)->expand();      const_cast<escript::Data*>(&in)->expand();
2065      paso::Coupler_startCollect(coupler, in.getSampleDataRO(0));      coupler->startCollect(in.getSampleDataRO(0));
2066    
2067      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2068      out.requireWrite();      out.requireWrite();
2069      const double* buffer = paso::Coupler_finishCollect(coupler);      const double* buffer = coupler->finishCollect();
2070    
2071  #pragma omp parallel for  #pragma omp parallel for
2072      for (index_t i=0; i<getNumNodes(); i++) {      for (index_t i=0; i<getNumNodes(); i++) {
# Line 2066  void Brick::dofToNodes(escript::Data& ou Line 2075  void Brick::dofToNodes(escript::Data& ou
2075                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
2076          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
2077      }      }
     paso::Coupler_free(coupler);  
2078  }  }
2079    
2080  //private  //private
# Line 2087  void Brick::populateSampleIds() Line 2095  void Brick::populateSampleIds()
2095          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2096      }      }
2097      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2098      m_nodeId.resize(getNumNodes());      
2099      m_dofId.resize(numDOF);      try {
2100      m_elementId.resize(getNumElements());          m_nodeId.resize(getNumNodes());
2101            m_dofId.resize(numDOF);
2102            m_elementId.resize(getNumElements());
2103        } catch (const std::length_error& le) {
2104            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2105        }
2106        
2107      // populate face element counts      // populate face element counts
2108      //left      //left
2109      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2569  void Brick::createPattern() Line 2582  void Brick::createPattern()
2582      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2583              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2584              &offsetInShared[0], 1, 0, m_mpiInfo));              &offsetInShared[0], 1, 0, m_mpiInfo));
2585      m_connector = paso::Connector_alloc(snd_shcomp, rcv_shcomp);      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
2586    
2587      // create main and couple blocks      // create main and couple blocks
2588      paso::Pattern *mainPattern = createMainPattern();      paso::Pattern_ptr mainPattern = createMainPattern();
2589      paso::Pattern *colPattern, *rowPattern;      paso::Pattern_ptr colPattern, rowPattern;
2590      createCouplePatterns(colIndices, rowIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, colPattern, rowPattern);
2591    
2592      // allocate paso distribution      // allocate paso distribution
2593      paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,      paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2594              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2595    
2596      // finally create the system matrix      // finally create the system matrix
2597      m_pattern = new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,      m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2598              distribution, distribution, mainPattern, colPattern, rowPattern,              distribution, distribution, mainPattern, colPattern, rowPattern,
2599              m_connector, m_connector);              m_connector, m_connector));
2600    
2601      // useful debug output      // useful debug output
2602      /*      /*
# Line 2642  void Brick::createPattern() Line 2655  void Brick::createPattern()
2655          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2656      }      }
2657      */      */
   
     paso::Pattern_free(mainPattern);  
     paso::Pattern_free(colPattern);  
     paso::Pattern_free(rowPattern);  
2658  }  }
2659    
2660  //private  //private
2661  void Brick::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,  void Brick::addToMatrixAndRHS(paso::SystemMatrix_ptr S, escript::Data& F,
2662           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2663           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2664  {  {
# Line 3333  int Brick::findNode(const double *coords Line 3342  int Brick::findNode(const double *coords
3342      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3343      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3344      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3345        
3346        //check if the point is even inside the domain
3347        if (x < 0 || y < 0 || z < 0
3348                || x > m_length[0] || y > m_length[1] || z > m_length[2])
3349            return NOT_MINE;
3350            
3351      // distance in elements      // distance in elements
3352      int ex = (int) floor(x / m_dx[0]);      int ex = (int) floor(x / m_dx[0]);
3353      int ey = (int) floor(y / m_dx[1]);      int ey = (int) floor(y / m_dx[1]);

Legend:
Removed from v.4816  
changed lines
  Added in v.4861

  ViewVC Help
Powered by ViewVC 1.1.26