/[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 4860 by sshaw, Thu Apr 10 04:08:42 2014 UTC
# Line 63  Brick::Brick(int n0, int n1, int n2, dou Line 63  Brick::Brick(int n0, int n1, int n2, dou
63               const simap_t& tagnamestonums) :               const simap_t& tagnamestonums) :
64      RipleyDomain(3)      RipleyDomain(3)
65  {  {
66        if (n0 <= 0 || n1 <= 0 || n2 <= 0)
67            throw RipleyException("Number of elements in each spatial dimension "
68                    "must be positive");
69    
70      // ignore subdivision parameters for serial run      // ignore subdivision parameters for serial run
71      if (m_mpiInfo->size == 1) {      if (m_mpiInfo->size == 1) {
72          d0=1;          d0=1;
# Line 203  Brick::Brick(int n0, int n1, int n2, dou Line 207  Brick::Brick(int n0, int n1, int n2, dou
207    
208  Brick::~Brick()  Brick::~Brick()
209  {  {
     paso::SystemMatrixPattern_free(m_pattern);  
     paso::Connector_free(m_connector);  
210      delete assembler;      delete assembler;
211  }  }
212    
# Line 2050  void Brick::nodesToDOF(escript::Data& ou Line 2052  void Brick::nodesToDOF(escript::Data& ou
2052  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const
2053  {  {
2054      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
2055      paso::Coupler* coupler = paso::Coupler_alloc(m_connector, numComp);      paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2056      // 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
2057      const_cast<escript::Data*>(&in)->expand();      const_cast<escript::Data*>(&in)->expand();
2058      paso::Coupler_startCollect(coupler, in.getSampleDataRO(0));      coupler->startCollect(in.getSampleDataRO(0));
2059    
2060      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2061      out.requireWrite();      out.requireWrite();
2062      const double* buffer = paso::Coupler_finishCollect(coupler);      const double* buffer = coupler->finishCollect();
2063    
2064  #pragma omp parallel for  #pragma omp parallel for
2065      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 2068  void Brick::dofToNodes(escript::Data& ou
2068                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
2069          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
2070      }      }
     paso::Coupler_free(coupler);  
2071  }  }
2072    
2073  //private  //private
# Line 2087  void Brick::populateSampleIds() Line 2088  void Brick::populateSampleIds()
2088          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2089      }      }
2090      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2091      m_nodeId.resize(getNumNodes());      
2092      m_dofId.resize(numDOF);      try {
2093      m_elementId.resize(getNumElements());          m_nodeId.resize(getNumNodes());
2094            m_dofId.resize(numDOF);
2095            m_elementId.resize(getNumElements());
2096        } catch (const std::length_error& le) {
2097            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2098        }
2099        
2100      // populate face element counts      // populate face element counts
2101      //left      //left
2102      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2569  void Brick::createPattern() Line 2575  void Brick::createPattern()
2575      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2576              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2577              &offsetInShared[0], 1, 0, m_mpiInfo));              &offsetInShared[0], 1, 0, m_mpiInfo));
2578      m_connector = paso::Connector_alloc(snd_shcomp, rcv_shcomp);      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
2579    
2580      // create main and couple blocks      // create main and couple blocks
2581      paso::Pattern *mainPattern = createMainPattern();      paso::Pattern_ptr mainPattern = createMainPattern();
2582      paso::Pattern *colPattern, *rowPattern;      paso::Pattern_ptr colPattern, rowPattern;
2583      createCouplePatterns(colIndices, rowIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, colPattern, rowPattern);
2584    
2585      // allocate paso distribution      // allocate paso distribution
2586      paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,      paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2587              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2588    
2589      // finally create the system matrix      // finally create the system matrix
2590      m_pattern = new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,      m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2591              distribution, distribution, mainPattern, colPattern, rowPattern,              distribution, distribution, mainPattern, colPattern, rowPattern,
2592              m_connector, m_connector);              m_connector, m_connector));
2593    
2594      // useful debug output      // useful debug output
2595      /*      /*
# Line 2642  void Brick::createPattern() Line 2648  void Brick::createPattern()
2648          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2649      }      }
2650      */      */
   
     paso::Pattern_free(mainPattern);  
     paso::Pattern_free(colPattern);  
     paso::Pattern_free(rowPattern);  
2651  }  }
2652    
2653  //private  //private
2654  void Brick::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,  void Brick::addToMatrixAndRHS(paso::SystemMatrix_ptr S, escript::Data& F,
2655           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2656           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2657  {  {
# Line 3333  int Brick::findNode(const double *coords Line 3335  int Brick::findNode(const double *coords
3335      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3336      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3337      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3338        
3339        //check if the point is even inside the domain
3340        if (x < 0 || y < 0 || z < 0
3341                || x > m_length[0] || y > m_length[1] || z > m_length[2])
3342            return NOT_MINE;
3343            
3344      // distance in elements      // distance in elements
3345      int ex = (int) floor(x / m_dx[0]);      int ex = (int) floor(x / m_dx[0]);
3346      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.4860

  ViewVC Help
Powered by ViewVC 1.1.26