/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

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

revision 3699 by caltinay, Thu Dec 1 22:59:42 2011 UTC revision 3707 by caltinay, Mon Dec 5 05:48:25 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    
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::interpolateOnDomain(escript::Data& target,  void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const
                                     const escript::Data& in) const  
258  {  {
259      const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));      escript::Data& in = *const_cast<escript::Data*>(&cIn);
260      const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));      const dim_t numComp = in.getDataPointSize();
261      if (inDomain != *this)      const double h0 = m_l0/m_gNE0;
262          throw RipleyException("Illegal domain of interpolant");      const double h1 = m_l1/m_gNE1;
263      if (targetDomain != *this)      if (out.getFunctionSpace().getTypeCode() == Elements) {
264          throw RipleyException("Illegal domain of interpolation target");          /* GENERATOR SNIP_GRAD_ELEMENTS TOP */
265            const double tmp0_13 = 0.78867513459481288225/h1;
266      stringstream msg;          const double tmp0_0 = 0.78867513459481288225/h0;
267      msg << "interpolateOnDomain() not implemented for function space "          const double tmp0_4 = 0.21132486540518711775/h0;
268          << in.getFunctionSpace().getTypeCode() << " -> "          const double tmp0_10 = 0.78867513459481288225/h1;
269          << target.getFunctionSpace().getTypeCode();          const double tmp0_1 = 0.21132486540518711775/h0;
270            const double tmp0_8 = -0.21132486540518711775/h1;
271      switch (in.getFunctionSpace().getTypeCode()) {          const double tmp0_14 = 0.21132486540518711775/h1;
272          case Nodes:          const double tmp0_5 = 0.78867513459481288225/h0;
273          case DegreesOfFreedom:          const double tmp0_11 = -0.78867513459481288225/h1;
274              switch (target.getFunctionSpace().getTypeCode()) {          const double tmp0_2 = -0.21132486540518711775/h0;
275                  case DegreesOfFreedom:          const double tmp0_9 = 0.21132486540518711775/h1;
276                      target=in;          const double tmp0_15 = -0.21132486540518711775/h1;
277                      break;          const double tmp0_6 = -0.78867513459481288225/h0;
278            const double tmp0_3 = -0.78867513459481288225/h0;
279                  case Elements:          const double tmp0_12 = -0.78867513459481288225/h1;
280                  {          const double tmp0_7 = -0.21132486540518711775/h0;
281                      const double tmp0_2 = 0.62200846792814621559;  #pragma omp parallel for
282                      const double tmp0_1 = 0.044658198738520451079;          for (index_t k1=0; k1 < m_NE1; ++k1) {
283                      const double tmp0_0 = 0.16666666666666666667;              for (index_t k0=0; k0 < m_NE0; ++k0) {
284                      const dim_t numComp = in.getDataPointSize();                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
285                      escript::Data* inn=const_cast<escript::Data*>(&in);                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
286  #pragma omp parallel for                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
287                      for (index_t k1=0; k1 < m_NE1; ++k1) {                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
288                          for (index_t k0=0; k0 < m_NE0; ++k0) {                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
289                              const register double* f_10 = inn->getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  for (index_t i=0; i < numComp; ++i) {
290                              const register double* f_11 = inn->getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                      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                              const register double* f_01 = inn->getSampleDataRO(INDEX2(k0,k1+1, m_N0));                      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                              const register double* f_00 = inn->getSampleDataRO(INDEX2(k0,k1, m_N0));                      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                              double* o = target.getSampleDataRW(INDEX2(k0,k1,m_NE0));                      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                              for (index_t i=0; i < numComp; ++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                                  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,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                                  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,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                                  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,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                                  o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);                  } /* end of component loop i */
299                              } // close component loop i              } /* end of k0 loop */
300                          } // close k0 loop          } /* end of k1 loop */
301                      } // close k1 loop          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
302                      break;      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
303                  }          /* GENERATOR SNIP_GRAD_FACES TOP */
304            if (m_faceOffset[0] > -1) {
305                  default:              const double tmp0_0 = 0.78867513459481288225/h0;
306                      throw RipleyException(msg.str());              const double tmp0_4 = 0.21132486540518711775/h0;
307              }              const double tmp0_1 = 0.21132486540518711775/h0;
308              break;              const double tmp0_8 = 1.0/h1;
309          default:              const double tmp0_5 = 0.78867513459481288225/h0;
310              throw RipleyException(msg.str());              const double tmp0_2 = -0.21132486540518711775/h0;
311                const double tmp0_9 = -1/h1;
312                const double tmp0_6 = -0.78867513459481288225/h0;
313                const double tmp0_3 = -0.78867513459481288225/h0;
314                const double tmp0_7 = -0.21132486540518711775/h0;
315    #pragma omp parallel for
316                for (index_t k1=0; k1 < m_NE1; ++k1) {
317                    const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
318                    const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
319                    const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
320                    const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
321                    double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
322                    for (index_t i=0; i < numComp; ++i) {
323                        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;
324                        o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;
325                        o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;
326                        o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;
327                    } /* end of component loop i */
328                } /* end of k1 loop */
329            } /* end of face 0 */
330            if (m_faceOffset[1] > -1) {
331                const double tmp0_0 = 0.78867513459481288225/h0;
332                const double tmp0_4 = 0.21132486540518711775/h0;
333                const double tmp0_1 = 0.21132486540518711775/h0;
334                const double tmp0_8 = -1/h1;
335                const double tmp0_5 = 0.78867513459481288225/h0;
336                const double tmp0_2 = -0.21132486540518711775/h0;
337                const double tmp0_9 = 1.0/h1;
338                const double tmp0_6 = -0.78867513459481288225/h0;
339                const double tmp0_3 = -0.78867513459481288225/h0;
340                const double tmp0_7 = -0.21132486540518711775/h0;
341    #pragma omp parallel for
342                for (index_t k1=0; k1 < m_NE1; ++k1) {
343                    const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
344                    const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
345                    const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
346                    const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
347                    double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
348                    for (index_t i=0; i < numComp; ++i) {
349                        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;
350                        o[INDEX3(i,1,0,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;
351                        o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;
352                        o[INDEX3(i,1,1,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;
353                    } /* end of component loop i */
354                } /* end of k1 loop */
355            } /* end of face 1 */
356            if (m_faceOffset[2] > -1) {
357                const double tmp0_0 = -1/h0;
358                const double tmp0_4 = 0.21132486540518711775/h1;
359                const double tmp0_1 = 1.0/h0;
360                const double tmp0_8 = 0.78867513459481288225/h1;
361                const double tmp0_5 = 0.78867513459481288225/h1;
362                const double tmp0_2 = -0.78867513459481288225/h1;
363                const double tmp0_9 = 0.21132486540518711775/h1;
364                const double tmp0_6 = -0.21132486540518711775/h1;
365                const double tmp0_3 = -0.21132486540518711775/h1;
366                const double tmp0_7 = -0.78867513459481288225/h1;
367    #pragma omp parallel for
368                for (index_t k0=0; k0 < m_NE0; ++k0) {
369                    const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
370                    const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
371                    const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
372                    const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
373                    double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
374                    for (index_t i=0; i < numComp; ++i) {
375                        o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
376                        o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_2 + f_01[i]*tmp0_5 + f_10[i]*tmp0_3 + f_11[i]*tmp0_4;
377                        o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
378                        o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_6 + f_01[i]*tmp0_9 + f_10[i]*tmp0_7 + f_11[i]*tmp0_8;
379                    } /* end of component loop i */
380                } /* end of k0 loop */
381            } /* end of face 2 */
382            if (m_faceOffset[3] > -1) {
383                const double tmp0_0 = 1.0/h0;
384                const double tmp0_4 = -0.21132486540518711775/h1;
385                const double tmp0_1 = -1/h0;
386                const double tmp0_8 = -0.78867513459481288225/h1;
387                const double tmp0_5 = -0.78867513459481288225/h1;
388                const double tmp0_2 = 0.21132486540518711775/h1;
389                const double tmp0_9 = -0.21132486540518711775/h1;
390                const double tmp0_6 = 0.78867513459481288225/h1;
391                const double tmp0_3 = 0.78867513459481288225/h1;
392                const double tmp0_7 = 0.21132486540518711775/h1;
393    #pragma omp parallel for
394                for (index_t k0=0; k0 < m_NE0; ++k0) {
395                    const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
396                    const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
397                    const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
398                    const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
399                    double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
400                    for (index_t i=0; i < numComp; ++i) {
401                        o[INDEX3(i,0,0,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
402                        o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_5 + f_01[i]*tmp0_3 + f_10[i]*tmp0_4 + f_11[i]*tmp0_2;
403                        o[INDEX3(i,0,1,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
404                        o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_7 + f_10[i]*tmp0_8 + f_11[i]*tmp0_6;
405                    } /* end of component loop i */
406                } /* end of k0 loop */
407            } /* end of face 3 */
408            /* GENERATOR SNIP_GRAD_FACES BOTTOM */
409        } else {
410            stringstream msg;
411            msg << "setToGradient() not implemented for "
412                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
413            throw RipleyException(msg.str());
414      }      }
415  }  }
416    
417    void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,
418            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
419            const escript::Data& C, const escript::Data& D,
420            const escript::Data& X, const escript::Data& Y,
421            const escript::Data& d, const escript::Data& y,
422            const escript::Data& d_contact, const escript::Data& y_contact,
423            const escript::Data& d_dirac, const escript::Data& y_dirac) const
424    {
425        throw RipleyException("addPDEToSystem() not implemented");
426    }
427    
428  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
429                                                  bool reducedColOrder) const                                                  bool reducedColOrder) const
430  {  {
# Line 397  Paso_SystemMatrixPattern* Rectangle::get Line 514  Paso_SystemMatrixPattern* Rectangle::get
514          recvShared.push_back(first);          recvShared.push_back(first);
515      }      }
516      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];
517        /*
518      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
519      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
520      for (size_t i=0; i<neighbour.size(); i++) {      for (size_t i=0; i<neighbour.size(); i++) {
# Line 410  Paso_SystemMatrixPattern* Rectangle::get Line 528  Paso_SystemMatrixPattern* Rectangle::get
528      for (size_t i=0; i<sendShared.size(); i++) {      for (size_t i=0; i<sendShared.size(); i++) {
529          cout << "shared[" << i << "]=" << sendShared[i] << endl;          cout << "shared[" << i << "]=" << sendShared[i] << endl;
530      }      }
531        */
532    
533      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
534              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
# Line 442  Paso_SystemMatrixPattern* Rectangle::get Line 561  Paso_SystemMatrixPattern* Rectangle::get
561      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
562              M, N, ptrC, indexC);              M, N, ptrC, indexC);
563    
564        /*
565      cout << "--- main_pattern ---" << endl;      cout << "--- main_pattern ---" << endl;
566      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << M << ", N=" << N << endl;
567      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<ptr.size(); i++) {
# Line 450  Paso_SystemMatrixPattern* Rectangle::get Line 570  Paso_SystemMatrixPattern* Rectangle::get
570      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<index.size(); i++) {
571          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << index[i] << endl;
572      }      }
573        */
574    
575      ptr.clear();      ptr.clear();
576      index.clear();      index.clear();
# Line 634  void Rectangle::populateSampleIds() Line 755  void Rectangle::populateSampleIds()
755          }          }
756      }      }
757    
758      // elements      // element id's
759      m_elementId.resize(getNumElements());      m_elementId.resize(getNumElements());
760  #pragma omp parallel for  #pragma omp parallel for
761      for (dim_t k=0; k<getNumElements(); k++) {      for (dim_t k=0; k<getNumElements(); k++) {
762          m_elementId[k]=k;          m_elementId[k]=k;
763      }      }
764    
765      // face elements      // face element id's
766      m_faceId.resize(getNumFaceElements());      m_faceId.resize(getNumFaceElements());
767  #pragma omp parallel for  #pragma omp parallel for
768      for (dim_t k=0; k<getNumFaceElements(); k++) {      for (dim_t k=0; k<getNumFaceElements(); k++) {
769          m_faceId[k]=k;          m_faceId[k]=k;
770      }      }
771    
772        // generate face offset vector
773        const IndexVector facesPerEdge = getNumFacesPerBoundary();
774        m_faceOffset.assign(facesPerEdge.size(), -1);
775        index_t offset=0;
776        for (size_t i=0; i<facesPerEdge.size(); i++) {
777            if (facesPerEdge[i]>0) {
778                m_faceOffset[i]=offset;
779                offset+=facesPerEdge[i];
780            }
781        }
782  }  }
783    
784  //private  //private
# Line 845  void Rectangle::generateCouplePatterns(P Line 977  void Rectangle::generateCouplePatterns(P
977          }          }
978      }      }
979    
980        /*
981      cout << "--- colCouple_pattern ---" << endl;      cout << "--- colCouple_pattern ---" << endl;
982      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << M << ", N=" << N << endl;
983      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<ptr.size(); i++) {
# Line 853  void Rectangle::generateCouplePatterns(P Line 986  void Rectangle::generateCouplePatterns(P
986      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<index.size(); i++) {
987          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << index[i] << endl;
988      }      }
989        */
990    
991      // now build the row couple pattern      // now build the row couple pattern
992      IndexVector ptr2(1,0);      IndexVector ptr2(1,0);
# Line 871  void Rectangle::generateCouplePatterns(P Line 1005  void Rectangle::generateCouplePatterns(P
1005          ptr2.push_back(ptr2.back()+n);          ptr2.push_back(ptr2.back()+n);
1006      }      }
1007    
1008        /*
1009      cout << "--- rowCouple_pattern ---" << endl;      cout << "--- rowCouple_pattern ---" << endl;
1010      cout << "M=" << N << ", N=" << M << endl;      cout << "M=" << N << ", N=" << M << endl;
1011      for (size_t i=0; i<ptr2.size(); i++) {      for (size_t i=0; i<ptr2.size(); i++) {
# Line 879  void Rectangle::generateCouplePatterns(P Line 1014  void Rectangle::generateCouplePatterns(P
1014      for (size_t i=0; i<index2.size(); i++) {      for (size_t i=0; i<index2.size(); i++) {
1015          cout << "index[" << i << "]=" << index2[i] << endl;          cout << "index[" << i << "]=" << index2[i] << endl;
1016      }      }
1017        */
1018    
1019      // paso will manage the memory      // paso will manage the memory
1020      index_t* indexC = MEMALLOC(index.size(), index_t);      index_t* indexC = MEMALLOC(index.size(), index_t);
# Line 895  void Rectangle::generateCouplePatterns(P Line 1031  void Rectangle::generateCouplePatterns(P
1031      *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);      *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);
1032  }  }
1033    
1034    //protected
1035    void Rectangle::interpolateNodesOnElements(escript::Data& out, escript::Data& in) const
1036    {
1037        const dim_t numComp = in.getDataPointSize();
1038        /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
1039        const double tmp0_2 = 0.62200846792814621559;
1040        const double tmp0_1 = 0.044658198738520451079;
1041        const double tmp0_0 = 0.16666666666666666667;
1042    #pragma omp parallel for
1043        for (index_t k1=0; k1 < m_NE1; ++k1) {
1044            for (index_t k0=0; k0 < m_NE0; ++k0) {
1045                const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1046                const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1047                const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1048                const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1049                double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1050                for (index_t i=0; i < numComp; ++i) {
1051                    o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);
1052                    o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);
1053                    o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);
1054                    o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);
1055                } /* end of component loop i */
1056            } /* end of k0 loop */
1057        } /* end of k1 loop */
1058        /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
1059    }
1060    
1061    //protected
1062    void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const
1063    {
1064        const dim_t numComp = in.getDataPointSize();
1065        /* GENERATOR SNIP_INTERPOLATE_FACES TOP */
1066        if (m_faceOffset[0] > -1) {
1067            const double tmp0_1 = 0.78867513459481288225;
1068            const double tmp0_0 = 0.21132486540518711775;
1069    #pragma omp parallel for
1070            for (index_t k1=0; k1 < m_NE1; ++k1) {
1071                const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1072                const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1073                double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1074                for (index_t i=0; i < numComp; ++i) {
1075                    o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_01[i]*tmp0_0;
1076                    o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_01[i]*tmp0_1;
1077                } /* end of component loop i */
1078            } /* end of k1 loop */
1079        } /* end of face 0 */
1080        if (m_faceOffset[1] > -1) {
1081            const double tmp0_1 = 0.21132486540518711775;
1082            const double tmp0_0 = 0.78867513459481288225;
1083    #pragma omp parallel for
1084            for (index_t k1=0; k1 < m_NE1; ++k1) {
1085                const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1086                const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1087                double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1088                for (index_t i=0; i < numComp; ++i) {
1089                    o[INDEX2(i,numComp,0)] = f_10[i]*tmp0_0 + f_11[i]*tmp0_1;
1090                    o[INDEX2(i,numComp,1)] = f_10[i]*tmp0_1 + f_11[i]*tmp0_0;
1091                } /* end of component loop i */
1092            } /* end of k1 loop */
1093        } /* end of face 1 */
1094        if (m_faceOffset[2] > -1) {
1095            const double tmp0_1 = 0.78867513459481288225;
1096            const double tmp0_0 = 0.21132486540518711775;
1097    #pragma omp parallel for
1098            for (index_t k0=0; k0 < m_NE0; ++k0) {
1099                const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1100                const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1101                double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1102                for (index_t i=0; i < numComp; ++i) {
1103                    o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_10[i]*tmp0_0;
1104                    o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
1105                } /* end of component loop i */
1106            } /* end of k0 loop */
1107        } /* end of face 2 */
1108        if (m_faceOffset[3] > -1) {
1109            const double tmp0_1 = 0.78867513459481288225;
1110            const double tmp0_0 = 0.21132486540518711775;
1111    #pragma omp parallel for
1112            for (index_t k0=0; k0 < m_NE0; ++k0) {
1113                const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1114                const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1115                double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1116                for (index_t i=0; i < numComp; ++i) {
1117                    o[INDEX2(i,numComp,0)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
1118                    o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_0 + f_11[i]*tmp0_1;
1119                } /* end of component loop i */
1120            } /* end of k0 loop */
1121        } /* end of face 3 */
1122        /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
1123    }
1124    
1125  } // end of namespace ripley  } // end of namespace ripley
1126    

Legend:
Removed from v.3699  
changed lines
  Added in v.3707

  ViewVC Help
Powered by ViewVC 1.1.26