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

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

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

revision 3703 by caltinay, Sun Dec 4 23:42:52 2011 UTC revision 3704 by caltinay, Mon Dec 5 01:59:08 2011 UTC
# Line 671  void Brick::populateSampleIds() Line 671  void Brick::populateSampleIds()
671      for (dim_t k=0; k<getNumFaceElements(); k++) {      for (dim_t k=0; k<getNumFaceElements(); k++) {
672          m_faceId[k]=k;          m_faceId[k]=k;
673      }      }
674    
675        // generate face offset vector
676        const IndexVector facesPerEdge = getNumFacesPerBoundary();
677        m_faceOffset.assign(facesPerEdge.size(), -1);
678        index_t offset=0;
679        for (size_t i=0; i<facesPerEdge.size(); i++) {
680            if (facesPerEdge[i]>0) {
681                m_faceOffset[i]=offset;
682                offset+=facesPerEdge[i];
683            }
684        }
685  }  }
686    
687  //protected  //protected
# Line 714  void Brick::interpolateNodesOnElements(e Line 725  void Brick::interpolateNodesOnElements(e
725  //protected  //protected
726  void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const  void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const
727  {  {
728      throw RipleyException("interpolateNodesOnFaces() not implemented");      const dim_t numComp = in.getDataPointSize();
729        /* GENERATOR SNIP_INTERPOLATE_FACES TOP */
730        if (m_faceOffset[0] > -1) {
731            const index_t k0 = 0;
732            const double tmp0_2 = 0.044658198738520451079;
733            const double tmp0_1 = 0.16666666666666666667;
734            const double tmp0_0 = 0.62200846792814621559;
735    #pragma omp parallel for
736            for (index_t k2=0; k2 < m_NE2; ++k2) {
737                for (index_t k1=0; k1 < m_NE1; ++k1) {
738                    const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
739                    const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
740                    const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
741                    const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
742                    double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
743                    for (index_t i=0; i < numComp; ++i) {
744                        o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_2 + tmp0_1*(f_001[i] + f_010[i]);
745                        o[INDEX2(i,numComp,1)] = f_001[i]*tmp0_2 + f_010[i]*tmp0_0 + tmp0_1*(f_000[i] + f_011[i]);
746                        o[INDEX2(i,numComp,2)] = f_001[i]*tmp0_0 + f_010[i]*tmp0_2 + tmp0_1*(f_000[i] + f_011[i]);
747                        o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_011[i]*tmp0_0 + tmp0_1*(f_001[i] + f_010[i]);
748                    } /* end of component loop i */
749                } /* end of k1 loop */
750            } /* end of k2 loop */
751        } /* end of face 0 */
752        if (m_faceOffset[1] > -1) {
753            const index_t k0 = 0;
754            const double tmp0_2 = 0.044658198738520451079;
755            const double tmp0_1 = 0.62200846792814621559;
756            const double tmp0_0 = 0.16666666666666666667;
757    #pragma omp parallel for
758            for (index_t k2=0; k2 < m_NE2; ++k2) {
759                for (index_t k1=0; k1 < m_NE1; ++k1) {
760                    const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
761                    const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
762                    const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
763                    const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
764                    double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
765                    for (index_t i=0; i < numComp; ++i) {
766                        o[INDEX2(i,numComp,0)] = f_100[i]*tmp0_1 + f_111[i]*tmp0_2 + tmp0_0*(f_101[i] + f_110[i]);
767                        o[INDEX2(i,numComp,1)] = f_101[i]*tmp0_2 + f_110[i]*tmp0_1 + tmp0_0*(f_100[i] + f_111[i]);
768                        o[INDEX2(i,numComp,2)] = f_101[i]*tmp0_1 + f_110[i]*tmp0_2 + tmp0_0*(f_100[i] + f_111[i]);
769                        o[INDEX2(i,numComp,3)] = f_100[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_101[i] + f_110[i]);
770                    } /* end of component loop i */
771                } /* end of k1 loop */
772            } /* end of k2 loop */
773        } /* end of face 1 */
774        if (m_faceOffset[2] > -1) {
775            const index_t k1 = 0;
776            const double tmp0_2 = 0.044658198738520451079;
777            const double tmp0_1 = 0.16666666666666666667;
778            const double tmp0_0 = 0.62200846792814621559;
779    #pragma omp parallel for
780            for (index_t k2=0; k2 < m_NE2; ++k2) {
781                for (index_t k0=0; k0 < m_NE0; ++k0) {
782                    const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
783                    const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
784                    const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
785                    const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
786                    double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
787                    for (index_t i=0; i < numComp; ++i) {
788                        o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_101[i]*tmp0_2 + tmp0_1*(f_001[i] + f_100[i]);
789                        o[INDEX2(i,numComp,1)] = f_001[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_101[i]);
790                        o[INDEX2(i,numComp,2)] = f_001[i]*tmp0_0 + f_100[i]*tmp0_2 + tmp0_1*(f_000[i] + f_101[i]);
791                        o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_100[i]);
792                    } /* end of component loop i */
793                } /* end of k0 loop */
794            } /* end of k2 loop */
795        } /* end of face 2 */
796        if (m_faceOffset[3] > -1) {
797            const index_t k1 = 0;
798            const double tmp0_2 = 0.044658198738520451079;
799            const double tmp0_1 = 0.62200846792814621559;
800            const double tmp0_0 = 0.16666666666666666667;
801    #pragma omp parallel for
802            for (index_t k2=0; k2 < m_NE2; ++k2) {
803                for (index_t k0=0; k0 < m_NE0; ++k0) {
804                    const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
805                    const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
806                    const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
807                    const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
808                    double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
809                    for (index_t i=0; i < numComp; ++i) {
810                        o[INDEX2(i,numComp,0)] = f_010[i]*tmp0_1 + f_111[i]*tmp0_2 + tmp0_0*(f_011[i] + f_110[i]);
811                        o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_2 + f_110[i]*tmp0_1 + tmp0_0*(f_010[i] + f_111[i]);
812                        o[INDEX2(i,numComp,2)] = f_011[i]*tmp0_1 + f_110[i]*tmp0_2 + tmp0_0*(f_010[i] + f_111[i]);
813                        o[INDEX2(i,numComp,3)] = f_010[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_011[i] + f_110[i]);
814                    } /* end of component loop i */
815                } /* end of k0 loop */
816            } /* end of k2 loop */
817        } /* end of face 3 */
818        if (m_faceOffset[4] > -1) {
819            const index_t k2 = 0;
820            const double tmp0_2 = 0.044658198738520451079;
821            const double tmp0_1 = 0.16666666666666666667;
822            const double tmp0_0 = 0.62200846792814621559;
823    #pragma omp parallel for
824            for (index_t k1=0; k1 < m_NE1; ++k1) {
825                for (index_t k0=0; k0 < m_NE0; ++k0) {
826                    const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
827                    const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
828                    const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
829                    const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
830                    double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
831                    for (index_t i=0; i < numComp; ++i) {
832                        o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_110[i]*tmp0_2 + tmp0_1*(f_010[i] + f_100[i]);
833                        o[INDEX2(i,numComp,1)] = f_010[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_110[i]);
834                        o[INDEX2(i,numComp,2)] = f_010[i]*tmp0_0 + f_100[i]*tmp0_2 + tmp0_1*(f_000[i] + f_110[i]);
835                        o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_110[i]*tmp0_0 + tmp0_1*(f_010[i] + f_100[i]);
836                    } /* end of component loop i */
837                } /* end of k0 loop */
838            } /* end of k1 loop */
839        } /* end of face 4 */
840        if (m_faceOffset[5] > -1) {
841            const index_t k2 = 0;
842            const double tmp0_2 = 0.044658198738520451079;
843            const double tmp0_1 = 0.16666666666666666667;
844            const double tmp0_0 = 0.62200846792814621559;
845    #pragma omp parallel for
846            for (index_t k1=0; k1 < m_NE1; ++k1) {
847                for (index_t k0=0; k0 < m_NE0; ++k0) {
848                    const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
849                    const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
850                    const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
851                    const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
852                    double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
853                    for (index_t i=0; i < numComp; ++i) {
854                        o[INDEX2(i,numComp,0)] = f_001[i]*tmp0_0 + f_111[i]*tmp0_2 + tmp0_1*(f_011[i] + f_101[i]);
855                        o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_111[i]);
856                        o[INDEX2(i,numComp,2)] = f_011[i]*tmp0_0 + f_101[i]*tmp0_2 + tmp0_1*(f_001[i] + f_111[i]);
857                        o[INDEX2(i,numComp,3)] = f_001[i]*tmp0_2 + f_111[i]*tmp0_0 + tmp0_1*(f_011[i] + f_101[i]);
858                    } /* end of component loop i */
859                } /* end of k0 loop */
860            } /* end of k1 loop */
861        } /* end of face 5 */
862        /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
863  }  }
864    
865  } // end of namespace ripley  } // end of namespace ripley

Legend:
Removed from v.3703  
changed lines
  Added in v.3704

  ViewVC Help
Powered by ViewVC 1.1.26