/[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 4941 by caltinay, Thu May 15 01:49:48 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 1999  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 2585  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 2659  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 2683  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.4940  
changed lines
  Added in v.4941

  ViewVC Help
Powered by ViewVC 1.1.26