/[escript]/branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3701 by caltinay, Fri Dec 2 02:04:58 2011 UTC revision 3702 by caltinay, Fri Dec 2 06:12:32 2011 UTC
# Line 239  const int* Rectangle::borrowSampleRefere Line 239  const int* Rectangle::borrowSampleRefere
239  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsCode, index_t id) const
240  {  {
241  #ifdef ESYS_MPI  #ifdef ESYS_MPI
242      if (fsCode == Nodes) {      if (fsCode != ReducedNodes) {
243          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
244          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
245          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
246      } else      } else {
247          throw RipleyException("ownSample() only implemented for Nodes");          stringstream msg;
248            msg << "ownSample() not implemented for "
249                << functionSpaceTypeAsString(fsCode);
250            throw RipleyException(msg.str());
251        }
252  #else  #else
253      return true;      return true;
254  #endif  #endif
255  }  }
256    
257  void Rectangle::interpolateNodesOnElements(escript::Data& out, escript::Data& in) const  void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const
258  {  {
259        escript::Data& in = *const_cast<escript::Data*>(&cIn);
260      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
261      /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */      const double h0 = m_l0/m_gNE0;
262      const double tmp0_2 = 0.62200846792814621559;      const double h1 = m_l1/m_gNE1;
263      const double tmp0_1 = 0.044658198738520451079;      if (out.getFunctionSpace().getTypeCode() == Elements) {
264      const double tmp0_0 = 0.16666666666666666667;          /* GENERATOR SNIP_GRAD_ELEMENTS TOP */
265            const double tmp0_13 = 0.78867513459481288225/h1;
266            const double tmp0_0 = 0.78867513459481288225/h0;
267            const double tmp0_4 = 0.21132486540518711775/h0;
268            const double tmp0_10 = 0.78867513459481288225/h1;
269            const double tmp0_1 = 0.21132486540518711775/h0;
270            const double tmp0_8 = -0.21132486540518711775/h1;
271            const double tmp0_14 = 0.21132486540518711775/h1;
272            const double tmp0_5 = 0.78867513459481288225/h0;
273            const double tmp0_11 = -0.78867513459481288225/h1;
274            const double tmp0_2 = -0.21132486540518711775/h0;
275            const double tmp0_9 = 0.21132486540518711775/h1;
276            const double tmp0_15 = -0.21132486540518711775/h1;
277            const double tmp0_6 = -0.78867513459481288225/h0;
278            const double tmp0_3 = -0.78867513459481288225/h0;
279            const double tmp0_12 = -0.78867513459481288225/h1;
280            const double tmp0_7 = -0.21132486540518711775/h0;
281  #pragma omp parallel for  #pragma omp parallel for
282      for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1 =0; k1 < m_NE1; ++k1) {
283          for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0 =0; k0 < m_NE0; ++k0) {
284              const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
285              const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
286              const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
287              const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
288              double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
289              for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
290                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2 + f_10[i]*tmp0_0 + f_11[i]*tmp0_1;
291                  o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_11 + f_01[i]*tmp0_10 + f_10[i]*tmp0_8 + f_11[i]*tmp0_9;
292                  o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2 + f_10[i]*tmp0_0 + f_11[i]*tmp0_1;
293                  o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_15 + f_01[i]*tmp0_14 + f_10[i]*tmp0_12 + f_11[i]*tmp0_13;
294              } /* end of component loop i */                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;
295          } /* end of k0 loop */                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*tmp0_11 + f_01[i]*tmp0_10 + f_10[i]*tmp0_8 + f_11[i]*tmp0_9;
296      } /* end of k1 loop */                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;
297      /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*tmp0_15 + f_01[i]*tmp0_14 + f_10[i]*tmp0_12 + f_11[i]*tmp0_13;
298  }                  } /* end of component loop i */
299                } /* end of k0 loop */
300  void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const          } /* end of k1 loop */
301  {          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
302      throw RipleyException("interpolateNodesOnFaces() not implemented");      } else {
303            throw RipleyException("setToGradient() not implemented");
304        }
305  }  }
306    
307  void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,  void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,
# Line 381  Paso_SystemMatrixPattern* Rectangle::get Line 404  Paso_SystemMatrixPattern* Rectangle::get
404          recvShared.push_back(first);          recvShared.push_back(first);
405      }      }
406      const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];      const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];
407        /*
408      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
409      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
410      for (size_t i=0; i<neighbour.size(); i++) {      for (size_t i=0; i<neighbour.size(); i++) {
# Line 394  Paso_SystemMatrixPattern* Rectangle::get Line 418  Paso_SystemMatrixPattern* Rectangle::get
418      for (size_t i=0; i<sendShared.size(); i++) {      for (size_t i=0; i<sendShared.size(); i++) {
419          cout << "shared[" << i << "]=" << sendShared[i] << endl;          cout << "shared[" << i << "]=" << sendShared[i] << endl;
420      }      }
421        */
422    
423      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
424              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
# Line 426  Paso_SystemMatrixPattern* Rectangle::get Line 451  Paso_SystemMatrixPattern* Rectangle::get
451      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
452              M, N, ptrC, indexC);              M, N, ptrC, indexC);
453    
454        /*
455      cout << "--- main_pattern ---" << endl;      cout << "--- main_pattern ---" << endl;
456      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << M << ", N=" << N << endl;
457      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<ptr.size(); i++) {
# Line 434  Paso_SystemMatrixPattern* Rectangle::get Line 460  Paso_SystemMatrixPattern* Rectangle::get
460      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<index.size(); i++) {
461          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << index[i] << endl;
462      }      }
463        */
464    
465      ptr.clear();      ptr.clear();
466      index.clear();      index.clear();
# Line 829  void Rectangle::generateCouplePatterns(P Line 856  void Rectangle::generateCouplePatterns(P
856          }          }
857      }      }
858    
859        /*
860      cout << "--- colCouple_pattern ---" << endl;      cout << "--- colCouple_pattern ---" << endl;
861      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << M << ", N=" << N << endl;
862      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<ptr.size(); i++) {
# Line 837  void Rectangle::generateCouplePatterns(P Line 865  void Rectangle::generateCouplePatterns(P
865      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<index.size(); i++) {
866          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << index[i] << endl;
867      }      }
868        */
869    
870      // now build the row couple pattern      // now build the row couple pattern
871      IndexVector ptr2(1,0);      IndexVector ptr2(1,0);
# Line 855  void Rectangle::generateCouplePatterns(P Line 884  void Rectangle::generateCouplePatterns(P
884          ptr2.push_back(ptr2.back()+n);          ptr2.push_back(ptr2.back()+n);
885      }      }
886    
887        /*
888      cout << "--- rowCouple_pattern ---" << endl;      cout << "--- rowCouple_pattern ---" << endl;
889      cout << "M=" << N << ", N=" << M << endl;      cout << "M=" << N << ", N=" << M << endl;
890      for (size_t i=0; i<ptr2.size(); i++) {      for (size_t i=0; i<ptr2.size(); i++) {
# Line 863  void Rectangle::generateCouplePatterns(P Line 893  void Rectangle::generateCouplePatterns(P
893      for (size_t i=0; i<index2.size(); i++) {      for (size_t i=0; i<index2.size(); i++) {
894          cout << "index[" << i << "]=" << index2[i] << endl;          cout << "index[" << i << "]=" << index2[i] << endl;
895      }      }
896        */
897    
898      // paso will manage the memory      // paso will manage the memory
899      index_t* indexC = MEMALLOC(index.size(), index_t);      index_t* indexC = MEMALLOC(index.size(), index_t);
# Line 879  void Rectangle::generateCouplePatterns(P Line 910  void Rectangle::generateCouplePatterns(P
910      *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);      *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);
911  }  }
912    
913    //protected
914    void Rectangle::interpolateNodesOnElements(escript::Data& out, escript::Data& in) const
915    {
916        const dim_t numComp = in.getDataPointSize();
917        /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
918        const double tmp0_2 = 0.62200846792814621559;
919        const double tmp0_1 = 0.044658198738520451079;
920        const double tmp0_0 = 0.16666666666666666667;
921    #pragma omp parallel for
922        for (index_t k1=0; k1 < m_NE1; ++k1) {
923            for (index_t k0=0; k0 < m_NE0; ++k0) {
924                const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
925                const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
926                const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
927                const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
928                double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
929                for (index_t i=0; i < numComp; ++i) {
930                    o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);
931                    o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);
932                    o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);
933                    o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);
934                } /* end of component loop i */
935            } /* end of k0 loop */
936        } /* end of k1 loop */
937        /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
938    }
939    
940    //protected
941    void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const
942    {
943        throw RipleyException("interpolateNodesOnFaces() not implemented");
944    }
945    
946  } // end of namespace ripley  } // end of namespace ripley
947    

Legend:
Removed from v.3701  
changed lines
  Added in v.3702

  ViewVC Help
Powered by ViewVC 1.1.26