/[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 3700 by caltinay, Thu Dec 1 22:59:42 2011 UTC revision 3701 by caltinay, Fri Dec 2 02:04:58 2011 UTC
# Line 232  const int* Rectangle::borrowSampleRefere Line 232  const int* Rectangle::borrowSampleRefere
232    
233      stringstream msg;      stringstream msg;
234      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs() not implemented for function space type "
235          << fsType;          << functionSpaceTypeAsString(fsType);
236      throw RipleyException(msg.str());      throw RipleyException(msg.str());
237  }  }
238    
# Line 250  bool Rectangle::ownSample(int fsCode, in Line 250  bool Rectangle::ownSample(int fsCode, in
250  #endif  #endif
251  }  }
252    
253  void Rectangle::interpolateOnDomain(escript::Data& target,  void Rectangle::interpolateNodesOnElements(escript::Data& out, escript::Data& in) const
                                     const escript::Data& in) const  
254  {  {
255      const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));      const dim_t numComp = in.getDataPointSize();
256      const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));      /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
257      if (inDomain != *this)      const double tmp0_2 = 0.62200846792814621559;
258          throw RipleyException("Illegal domain of interpolant");      const double tmp0_1 = 0.044658198738520451079;
259      if (targetDomain != *this)      const double tmp0_0 = 0.16666666666666666667;
         throw RipleyException("Illegal domain of interpolation target");  
   
     stringstream msg;  
     msg << "interpolateOnDomain() not implemented for function space "  
         << in.getFunctionSpace().getTypeCode() << " -> "  
         << target.getFunctionSpace().getTypeCode();  
   
     switch (in.getFunctionSpace().getTypeCode()) {  
         case Nodes:  
         case DegreesOfFreedom:  
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case DegreesOfFreedom:  
                     target=in;  
                     break;  
   
                 case Elements:  
                 {  
                     const double tmp0_2 = 0.62200846792814621559;  
                     const double tmp0_1 = 0.044658198738520451079;  
                     const double tmp0_0 = 0.16666666666666666667;  
                     const dim_t numComp = in.getDataPointSize();  
                     escript::Data* inn=const_cast<escript::Data*>(&in);  
260  #pragma omp parallel for  #pragma omp parallel for
261                      for (index_t k1=0; k1 < m_NE1; ++k1) {      for (index_t k1=0; k1 < m_NE1; ++k1) {
262                          for (index_t k0=0; k0 < m_NE0; ++k0) {          for (index_t k0=0; k0 < m_NE0; ++k0) {
263                              const register double* f_10 = inn->getSampleDataRO(INDEX2(k0+1,k1, m_N0));              const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
264                              const register double* f_11 = inn->getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
265                              const register double* f_01 = inn->getSampleDataRO(INDEX2(k0,k1+1, m_N0));              const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
266                              const register double* f_00 = inn->getSampleDataRO(INDEX2(k0,k1, m_N0));              const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
267                              double* o = target.getSampleDataRW(INDEX2(k0,k1,m_NE0));              double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
268                              for (index_t i=0; i < numComp; ++i) {              for (index_t i=0; i < numComp; ++i) {
269                                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);
270                                  o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);                  o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);
271                                  o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);                  o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);
272                                  o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);                  o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);
273                              } // close component loop i              } /* end of component loop i */
274                          } // close k0 loop          } /* end of k0 loop */
275                      } // close k1 loop      } /* end of k1 loop */
276                      break;      /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
277                  }  }
278    
279                  default:  void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const
280                      throw RipleyException(msg.str());  {
281              }      throw RipleyException("interpolateNodesOnFaces() not implemented");
282              break;  }
283          default:  
284              throw RipleyException(msg.str());  void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,
285      }          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
286            const escript::Data& C, const escript::Data& D,
287            const escript::Data& X, const escript::Data& Y,
288            const escript::Data& d, const escript::Data& y,
289            const escript::Data& d_contact, const escript::Data& y_contact,
290            const escript::Data& d_dirac, const escript::Data& y_dirac) const
291    {
292        throw RipleyException("addPDEToSystem() not implemented");
293  }  }
294    
295  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
Legend:
Removed from v.3700  
changed lines
  Added in v.3701
  ViewVC Help
Powered by ViewVC 1.1.26