/[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 3702 by caltinay, Fri Dec 2 06:12:32 2011 UTC revision 3704 by caltinay, Mon Dec 5 01:59:08 2011 UTC
# Line 645  void Rectangle::populateSampleIds() Line 645  void Rectangle::populateSampleIds()
645          }          }
646      }      }
647    
648      // elements      // element id's
649      m_elementId.resize(getNumElements());      m_elementId.resize(getNumElements());
650  #pragma omp parallel for  #pragma omp parallel for
651      for (dim_t k=0; k<getNumElements(); k++) {      for (dim_t k=0; k<getNumElements(); k++) {
652          m_elementId[k]=k;          m_elementId[k]=k;
653      }      }
654    
655      // face elements      // face element id's
656      m_faceId.resize(getNumFaceElements());      m_faceId.resize(getNumFaceElements());
657  #pragma omp parallel for  #pragma omp parallel for
658      for (dim_t k=0; k<getNumFaceElements(); k++) {      for (dim_t k=0; k<getNumFaceElements(); k++) {
659          m_faceId[k]=k;          m_faceId[k]=k;
660      }      }
661    
662        // generate face offset vector
663        const IndexVector facesPerEdge = getNumFacesPerBoundary();
664        m_faceOffset.assign(facesPerEdge.size(), -1);
665        index_t offset=0;
666        for (size_t i=0; i<facesPerEdge.size(); i++) {
667            if (facesPerEdge[i]>0) {
668                m_faceOffset[i]=offset;
669                offset+=facesPerEdge[i];
670            }
671        }
672  }  }
673    
674  //private  //private
# Line 940  void Rectangle::interpolateNodesOnElemen Line 951  void Rectangle::interpolateNodesOnElemen
951  //protected  //protected
952  void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const  void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const
953  {  {
954      throw RipleyException("interpolateNodesOnFaces() not implemented");      const dim_t numComp = in.getDataPointSize();
955        /* GENERATOR SNIP_INTERPOLATE_FACES TOP */
956        if (m_faceOffset[0] > -1) {
957            const index_t k0 = 0;
958            const double tmp0_1 = 0.78867513459481288225;
959            const double tmp0_0 = 0.21132486540518711775;
960    #pragma omp parallel for
961            for (index_t k1=0; k1 < m_NE1; ++k1) {
962                const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
963                const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
964                double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k0,k1,m_NE0));
965                for (index_t i=0; i < numComp; ++i) {
966                    o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_01[i]*tmp0_0;
967                    o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_01[i]*tmp0_1;
968                } /* end of component loop i */
969            } /* end of k1 loop */
970        } /* end of face 0 */
971        if (m_faceOffset[1] > -1) {
972            const index_t k0 = 0;
973            const double tmp0_1 = 0.21132486540518711775;
974            const double tmp0_0 = 0.78867513459481288225;
975    #pragma omp parallel for
976            for (index_t k1=0; k1 < m_NE1; ++k1) {
977                const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
978                const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
979                double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k0,k1,m_NE0));
980                for (index_t i=0; i < numComp; ++i) {
981                    o[INDEX2(i,numComp,0)] = f_10[i]*tmp0_0 + f_11[i]*tmp0_1;
982                    o[INDEX2(i,numComp,1)] = f_10[i]*tmp0_1 + f_11[i]*tmp0_0;
983                } /* end of component loop i */
984            } /* end of k1 loop */
985        } /* end of face 1 */
986        if (m_faceOffset[2] > -1) {
987            const index_t k1 = 0;
988            const double tmp0_1 = 0.78867513459481288225;
989            const double tmp0_0 = 0.21132486540518711775;
990    #pragma omp parallel for
991            for (index_t k0=0; k0 < m_NE0; ++k0) {
992                const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
993                const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
994                double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k1,m_NE0));
995                for (index_t i=0; i < numComp; ++i) {
996                    o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_10[i]*tmp0_0;
997                    o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
998                } /* end of component loop i */
999            } /* end of k0 loop */
1000        } /* end of face 2 */
1001        if (m_faceOffset[3] > -1) {
1002            const index_t k1 = 0;
1003            const double tmp0_1 = 0.78867513459481288225;
1004            const double tmp0_0 = 0.21132486540518711775;
1005    #pragma omp parallel for
1006            for (index_t k0=0; k0 < m_NE0; ++k0) {
1007                const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1008                const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1009                double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k1,m_NE0));
1010                for (index_t i=0; i < numComp; ++i) {
1011                    o[INDEX2(i,numComp,0)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
1012                    o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_0 + f_11[i]*tmp0_1;
1013                } /* end of component loop i */
1014            } /* end of k0 loop */
1015        } /* end of face 3 */
1016        /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
1017  }  }
1018    
1019  } // end of namespace ripley  } // end of namespace ripley
Legend:
Removed from v.3702  
changed lines
  Added in v.3704
  ViewVC Help
Powered by ViewVC 1.1.26