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

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

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

revision 3710 by caltinay, Mon Dec 5 05:48:25 2011 UTC revision 3711 by caltinay, Tue Dec 6 00:24:43 2011 UTC
# Line 373  void Brick::setToGradient(escript::Data& Line 373  void Brick::setToGradient(escript::Data&
373              } /* end of k1 loop */              } /* end of k1 loop */
374          } /* end of k2 loop */          } /* end of k2 loop */
375          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
376        } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
377            /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */
378            const double tmp0_0 = -0.25/h0;
379            const double tmp0_4 = -0.25/h2;
380            const double tmp0_1 = 0.25/h0;
381            const double tmp0_5 = 0.25/h2;
382            const double tmp0_2 = -0.25/h1;
383            const double tmp0_3 = 0.25/h1;
384    #pragma omp parallel for
385            for (index_t k2=0; k2 < m_NE2; ++k2) {
386                for (index_t k1=0; k1 < m_NE1; ++k1) {
387                    for (index_t k0=0; k0 < m_NE0; ++k0) {
388                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
389                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
390                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
391                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
392                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
393                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
394                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
395                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
396                        double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
397                        for (index_t i=0; i < numComp; ++i) {
398                            o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]) + tmp0_1*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
399                            o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_001[i] + f_100[i] + f_101[i]) + tmp0_3*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);
400                            o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_010[i] + f_100[i] + f_110[i]) + tmp0_5*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);
401                        } /* end of component loop i */
402                    } /* end of k0 loop */
403                } /* end of k1 loop */
404            } /* end of k2 loop */
405            /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
406      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
407          /* GENERATOR SNIP_GRAD_FACES TOP */          /* GENERATOR SNIP_GRAD_FACES TOP */
408          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
# Line 772  void Brick::setToGradient(escript::Data& Line 802  void Brick::setToGradient(escript::Data&
802              } /* end of k1 loop */              } /* end of k1 loop */
803          } /* end of face 5 */          } /* end of face 5 */
804          /* GENERATOR SNIP_GRAD_FACES BOTTOM */          /* GENERATOR SNIP_GRAD_FACES BOTTOM */
805        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
806            /* GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
807            if (m_faceOffset[0] > -1) {
808                const double tmp0_0 = -0.25/h0;
809                const double tmp0_4 = -0.5/h2;
810                const double tmp0_1 = 0.25/h0;
811                const double tmp0_5 = 0.5/h2;
812                const double tmp0_2 = -0.5/h1;
813                const double tmp0_3 = 0.5/h1;
814    #pragma omp parallel for
815                for (index_t k2=0; k2 < m_NE2; ++k2) {
816                    for (index_t k1=0; k1 < m_NE1; ++k1) {
817                        const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
818                        const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
819                        const register double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));
820                        const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
821                        const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
822                        const register double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));
823                        const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
824                        const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
825                        double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
826                        for (index_t i=0; i < numComp; ++i) {
827                            o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]) + tmp0_1*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
828                            o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_001[i]) + tmp0_3*(f_010[i] + f_011[i]);
829                            o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_010[i]) + tmp0_5*(f_001[i] + f_011[i]);
830                        } /* end of component loop i */
831                    } /* end of k1 loop */
832                } /* end of k2 loop */
833            } /* end of face 0 */
834            if (m_faceOffset[1] > -1) {
835                const double tmp0_0 = -0.25/h0;
836                const double tmp0_4 = 0.5/h2;
837                const double tmp0_1 = 0.25/h0;
838                const double tmp0_5 = -0.5/h2;
839                const double tmp0_2 = -0.5/h1;
840                const double tmp0_3 = 0.5/h1;
841    #pragma omp parallel for
842                for (index_t k2=0; k2 < m_NE2; ++k2) {
843                    for (index_t k1=0; k1 < m_NE1; ++k1) {
844                        const register double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));
845                        const register double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));
846                        const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
847                        const register double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));
848                        const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
849                        const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
850                        const register double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));
851                        const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
852                        double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
853                        for (index_t i=0; i < numComp; ++i) {
854                            o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]) + tmp0_1*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
855                            o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_100[i] + f_101[i]) + tmp0_3*(f_110[i] + f_111[i]);
856                            o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_101[i] + f_111[i]) + tmp0_5*(f_100[i] + f_110[i]);
857                        } /* end of component loop i */
858                    } /* end of k1 loop */
859                } /* end of k2 loop */
860            } /* end of face 1 */
861            if (m_faceOffset[2] > -1) {
862                const double tmp0_0 = -0.5/h0;
863                const double tmp0_4 = -0.5/h2;
864                const double tmp0_1 = 0.5/h0;
865                const double tmp0_5 = 0.5/h2;
866                const double tmp0_2 = -0.25/h1;
867                const double tmp0_3 = 0.25/h1;
868    #pragma omp parallel for
869                for (index_t k2=0; k2 < m_NE2; ++k2) {
870                    for (index_t k0=0; k0 < m_NE0; ++k0) {
871                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
872                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
873                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
874                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
875                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));
876                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));
877                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
878                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
879                        double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
880                        for (index_t i=0; i < numComp; ++i) {
881                            o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_001[i]) + tmp0_1*(f_100[i] + f_101[i]);
882                            o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_001[i] + f_100[i] + f_101[i]) + tmp0_3*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);
883                            o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_100[i]) + tmp0_5*(f_001[i] + f_101[i]);
884                        } /* end of component loop i */
885                    } /* end of k0 loop */
886                } /* end of k2 loop */
887            } /* end of face 2 */
888            if (m_faceOffset[3] > -1) {
889                const double tmp0_0 = -0.5/h0;
890                const double tmp0_4 = 0.5/h2;
891                const double tmp0_1 = 0.5/h0;
892                const double tmp0_5 = -0.5/h2;
893                const double tmp0_2 = 0.25/h1;
894                const double tmp0_3 = -0.25/h1;
895    #pragma omp parallel for
896                for (index_t k2=0; k2 < m_NE2; ++k2) {
897                    for (index_t k0=0; k0 < m_NE0; ++k0) {
898                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
899                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
900                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
901                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
902                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));
903                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));
904                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));
905                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
906                        double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
907                        for (index_t i=0; i < numComp; ++i) {
908                            o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_010[i] + f_011[i]) + tmp0_1*(f_110[i] + f_111[i]);
909                            o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_010[i] + f_011[i] + f_110[i] + f_111[i]) + tmp0_3*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);
910                            o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_011[i] + f_111[i]) + tmp0_5*(f_010[i] + f_110[i]);
911                        } /* end of component loop i */
912                    } /* end of k0 loop */
913                } /* end of k2 loop */
914            } /* end of face 3 */
915            if (m_faceOffset[4] > -1) {
916                const double tmp0_0 = -0.5/h0;
917                const double tmp0_4 = -0.25/h2;
918                const double tmp0_1 = 0.5/h0;
919                const double tmp0_5 = 0.25/h2;
920                const double tmp0_2 = -0.5/h1;
921                const double tmp0_3 = 0.5/h1;
922    #pragma omp parallel for
923                for (index_t k1=0; k1 < m_NE1; ++k1) {
924                    for (index_t k0=0; k0 < m_NE0; ++k0) {
925                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
926                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
927                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
928                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
929                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));
930                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));
931                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));
932                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
933                        double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
934                        for (index_t i=0; i < numComp; ++i) {
935                            o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_010[i]) + tmp0_1*(f_100[i] + f_110[i]);
936                            o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_100[i]) + tmp0_3*(f_010[i] + f_110[i]);
937                            o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_010[i] + f_100[i] + f_110[i]) + tmp0_5*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);
938                        } /* end of component loop i */
939                    } /* end of k0 loop */
940                } /* end of k1 loop */
941            } /* end of face 4 */
942            if (m_faceOffset[5] > -1) {
943                const double tmp0_0 = -0.5/h0;
944                const double tmp0_4 = 0.25/h2;
945                const double tmp0_1 = 0.5/h0;
946                const double tmp0_5 = -0.25/h2;
947                const double tmp0_2 = -0.5/h1;
948                const double tmp0_3 = 0.5/h1;
949    #pragma omp parallel for
950                for (index_t k1=0; k1 < m_NE1; ++k1) {
951                    for (index_t k0=0; k0 < m_NE0; ++k0) {
952                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
953                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
954                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
955                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
956                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));
957                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
958                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));
959                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
960                        double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
961                        for (index_t i=0; i < numComp; ++i) {
962                            o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_001[i] + f_011[i]) + tmp0_1*(f_101[i] + f_111[i]);
963                            o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_001[i] + f_101[i]) + tmp0_3*(f_011[i] + f_111[i]);
964                            o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_001[i] + f_011[i] + f_101[i] + f_111[i]) + tmp0_5*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);
965                        } /* end of component loop i */
966                    } /* end of k0 loop */
967                } /* end of k1 loop */
968            } /* end of face 5 */
969            /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */
970      } else {      } else {
971          stringstream msg;          stringstream msg;
972          msg << "setToGradient() not implemented for "          msg << "setToGradient() not implemented for "
# Line 1087  void Brick::populateSampleIds() Line 1282  void Brick::populateSampleIds()
1282  }  }
1283    
1284  //protected  //protected
1285  void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in) const  void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in,
1286                                           bool reduced) const
1287  {  {
1288      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1289      /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */      if (reduced) {
1290      const double tmp0_3 = 0.0094373878376559314545;          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */
1291      const double tmp0_2 = 0.035220810900864519624;          const double tmp0_0 = 0.12500000000000000000;
1292      const double tmp0_1 = 0.13144585576580214704;  #pragma omp parallel for
1293      const double tmp0_0 = 0.49056261216234406855;          for (index_t k2=0; k2 < m_NE2; ++k2) {
1294  #pragma omp parallel for              for (index_t k1=0; k1 < m_NE1; ++k1) {
1295      for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1296          for (index_t k1=0; k1 < m_NE1; ++k1) {                      const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
1297              for (index_t k0=0; k0 < m_NE0; ++k0) {                      const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
1298                  const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));                      const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
1299                  const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));                      const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
1300                  const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));                      const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
1301                  const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));                      const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
1302                  const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));                      const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
1303                  const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));                      const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1304                  const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1305                  const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));                      for (index_t i=0; i < numComp; ++i) {
1306                  double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                          o[INDEX2(i,numComp,0)] = tmp0_0*(f_000[i] + f_001[i] + f_010[i] + f_011[i] + f_100[i] + f_101[i] + f_110[i] + f_111[i]);
1307                  for (index_t i=0; i < numComp; ++i) {                      } /* end of component loop i */
1308                      o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i] + f_100[i]) + tmp0_2*(f_011[i] + f_101[i] + f_110[i]);                  } /* end of k0 loop */
1309                      o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_3 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_101[i] + f_110[i]) + tmp0_2*(f_001[i] + f_010[i] + f_111[i]);              } /* end of k1 loop */
1310                      o[INDEX2(i,numComp,2)] = f_010[i]*tmp0_0 + f_101[i]*tmp0_3 + tmp0_1*(f_000[i] + f_011[i] + f_110[i]) + tmp0_2*(f_001[i] + f_100[i] + f_111[i]);          } /* end of k2 loop */
1311                      o[INDEX2(i,numComp,3)] = f_001[i]*tmp0_3 + f_110[i]*tmp0_0 + tmp0_1*(f_010[i] + f_100[i] + f_111[i]) + tmp0_2*(f_000[i] + f_011[i] + f_101[i]);          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */
1312                      o[INDEX2(i,numComp,4)] = f_001[i]*tmp0_0 + f_110[i]*tmp0_3 + tmp0_1*(f_000[i] + f_011[i] + f_101[i]) + tmp0_2*(f_010[i] + f_100[i] + f_111[i]);      } else {
1313                      o[INDEX2(i,numComp,5)] = f_010[i]*tmp0_3 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_100[i] + f_111[i]) + tmp0_2*(f_000[i] + f_011[i] + f_110[i]);          /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
1314                      o[INDEX2(i,numComp,6)] = f_011[i]*tmp0_0 + f_100[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i] + f_111[i]) + tmp0_2*(f_000[i] + f_101[i] + f_110[i]);          const double tmp0_3 = 0.0094373878376559314545;
1315                      o[INDEX2(i,numComp,7)] = f_000[i]*tmp0_3 + f_111[i]*tmp0_0 + tmp0_1*(f_011[i] + f_101[i] + f_110[i]) + tmp0_2*(f_001[i] + f_010[i] + f_100[i]);          const double tmp0_2 = 0.035220810900864519624;
1316                  } /* end of component loop i */          const double tmp0_1 = 0.13144585576580214704;
1317              } /* end of k0 loop */          const double tmp0_0 = 0.49056261216234406855;
1318          } /* end of k1 loop */  #pragma omp parallel for
1319      } /* end of k2 loop */          for (index_t k2=0; k2 < m_NE2; ++k2) {
1320      /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */              for (index_t k1=0; k1 < m_NE1; ++k1) {
1321                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1322                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
1323                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
1324                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
1325                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
1326                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
1327                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
1328                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
1329                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1330                        double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1331                        for (index_t i=0; i < numComp; ++i) {
1332                            o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i] + f_100[i]) + tmp0_2*(f_011[i] + f_101[i] + f_110[i]);
1333                            o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_3 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_101[i] + f_110[i]) + tmp0_2*(f_001[i] + f_010[i] + f_111[i]);
1334                            o[INDEX2(i,numComp,2)] = f_010[i]*tmp0_0 + f_101[i]*tmp0_3 + tmp0_1*(f_000[i] + f_011[i] + f_110[i]) + tmp0_2*(f_001[i] + f_100[i] + f_111[i]);
1335                            o[INDEX2(i,numComp,3)] = f_001[i]*tmp0_3 + f_110[i]*tmp0_0 + tmp0_1*(f_010[i] + f_100[i] + f_111[i]) + tmp0_2*(f_000[i] + f_011[i] + f_101[i]);
1336                            o[INDEX2(i,numComp,4)] = f_001[i]*tmp0_0 + f_110[i]*tmp0_3 + tmp0_1*(f_000[i] + f_011[i] + f_101[i]) + tmp0_2*(f_010[i] + f_100[i] + f_111[i]);
1337                            o[INDEX2(i,numComp,5)] = f_010[i]*tmp0_3 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_100[i] + f_111[i]) + tmp0_2*(f_000[i] + f_011[i] + f_110[i]);
1338                            o[INDEX2(i,numComp,6)] = f_011[i]*tmp0_0 + f_100[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i] + f_111[i]) + tmp0_2*(f_000[i] + f_101[i] + f_110[i]);
1339                            o[INDEX2(i,numComp,7)] = f_000[i]*tmp0_3 + f_111[i]*tmp0_0 + tmp0_1*(f_011[i] + f_101[i] + f_110[i]) + tmp0_2*(f_001[i] + f_010[i] + f_100[i]);
1340                        } /* end of component loop i */
1341                    } /* end of k0 loop */
1342                } /* end of k1 loop */
1343            } /* end of k2 loop */
1344            /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
1345        }
1346  }  }
1347    
1348  //protected  //protected
1349  void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const  void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,
1350                                        bool reduced) const
1351  {  {
1352      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1353      /* GENERATOR SNIP_INTERPOLATE_FACES TOP */      if (reduced) {
1354      if (m_faceOffset[0] > -1) {          /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */
1355          const double tmp0_2 = 0.044658198738520451079;          if (m_faceOffset[0] > -1) {
1356          const double tmp0_1 = 0.16666666666666666667;              const double tmp0_0 = 0.25000000000000000000;
1357          const double tmp0_0 = 0.62200846792814621559;  #pragma omp parallel for
1358                for (index_t k2=0; k2 < m_NE2; ++k2) {
1359                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1360                        const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
1361                        const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
1362                        const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
1363                        const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
1364                        double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1365                        for (index_t i=0; i < numComp; ++i) {
1366                            o[INDEX2(i,numComp,0)] = tmp0_0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]);
1367                        } /* end of component loop i */
1368                    } /* end of k1 loop */
1369                } /* end of k2 loop */
1370            } /* end of face 0 */
1371            if (m_faceOffset[1] > -1) {
1372                const double tmp0_0 = 0.25000000000000000000;
1373    #pragma omp parallel for
1374                for (index_t k2=0; k2 < m_NE2; ++k2) {
1375                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1376                        const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
1377                        const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
1378                        const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
1379                        const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
1380                        double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1381                        for (index_t i=0; i < numComp; ++i) {
1382                            o[INDEX2(i,numComp,0)] = tmp0_0*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
1383                        } /* end of component loop i */
1384                    } /* end of k1 loop */
1385                } /* end of k2 loop */
1386            } /* end of face 1 */
1387            if (m_faceOffset[2] > -1) {
1388                const double tmp0_0 = 0.25000000000000000000;
1389    #pragma omp parallel for
1390                for (index_t k2=0; k2 < m_NE2; ++k2) {
1391                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1392                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
1393                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
1394                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
1395                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
1396                        double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1397                        for (index_t i=0; i < numComp; ++i) {
1398                            o[INDEX2(i,numComp,0)] = tmp0_0*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);
1399                        } /* end of component loop i */
1400                    } /* end of k0 loop */
1401                } /* end of k2 loop */
1402            } /* end of face 2 */
1403            if (m_faceOffset[3] > -1) {
1404                const double tmp0_0 = 0.25000000000000000000;
1405    #pragma omp parallel for
1406                for (index_t k2=0; k2 < m_NE2; ++k2) {
1407                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1408                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
1409                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
1410                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
1411                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
1412                        double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1413                        for (index_t i=0; i < numComp; ++i) {
1414                            o[INDEX2(i,numComp,0)] = tmp0_0*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);
1415                        } /* end of component loop i */
1416                    } /* end of k0 loop */
1417                } /* end of k2 loop */
1418            } /* end of face 3 */
1419            if (m_faceOffset[4] > -1) {
1420                const double tmp0_0 = 0.25000000000000000000;
1421  #pragma omp parallel for  #pragma omp parallel for
         for (index_t k2=0; k2 < m_NE2; ++k2) {  
1422              for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k1=0; k1 < m_NE1; ++k1) {
1423                  const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1424                  const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));                      const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
1425                  const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));                      const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
1426                  const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));                      const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
1427                  double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                      const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
1428                  for (index_t i=0; i < numComp; ++i) {                      double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1429                      o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_2 + tmp0_1*(f_001[i] + f_010[i]);                      for (index_t i=0; i < numComp; ++i) {
1430                      o[INDEX2(i,numComp,1)] = f_001[i]*tmp0_2 + f_010[i]*tmp0_0 + tmp0_1*(f_000[i] + f_011[i]);                          o[INDEX2(i,numComp,0)] = tmp0_0*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);
1431                      o[INDEX2(i,numComp,2)] = f_001[i]*tmp0_0 + f_010[i]*tmp0_2 + tmp0_1*(f_000[i] + f_011[i]);                      } /* end of component loop i */
1432                      o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_011[i]*tmp0_0 + tmp0_1*(f_001[i] + f_010[i]);                  } /* end of k0 loop */
                 } /* end of component loop i */  
1433              } /* end of k1 loop */              } /* end of k1 loop */
1434          } /* end of k2 loop */          } /* end of face 4 */
1435      } /* end of face 0 */          if (m_faceOffset[5] > -1) {
1436      if (m_faceOffset[1] > -1) {              const double tmp0_0 = 0.25000000000000000000;
         const double tmp0_2 = 0.044658198738520451079;  
         const double tmp0_1 = 0.62200846792814621559;  
         const double tmp0_0 = 0.16666666666666666667;  
1437  #pragma omp parallel for  #pragma omp parallel for
         for (index_t k2=0; k2 < m_NE2; ++k2) {  
1438              for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k1=0; k1 < m_NE1; ++k1) {
1439                  const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1440                  const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));                      const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
1441                  const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));                      const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
1442                  const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));                      const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
1443                  double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                      const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
1444                  for (index_t i=0; i < numComp; ++i) {                      double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1445                      o[INDEX2(i,numComp,0)] = f_100[i]*tmp0_1 + f_111[i]*tmp0_2 + tmp0_0*(f_101[i] + f_110[i]);                      for (index_t i=0; i < numComp; ++i) {
1446                      o[INDEX2(i,numComp,1)] = f_101[i]*tmp0_2 + f_110[i]*tmp0_1 + tmp0_0*(f_100[i] + f_111[i]);                          o[INDEX2(i,numComp,0)] = tmp0_0*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);
1447                      o[INDEX2(i,numComp,2)] = f_101[i]*tmp0_1 + f_110[i]*tmp0_2 + tmp0_0*(f_100[i] + f_111[i]);                      } /* end of component loop i */
1448                      o[INDEX2(i,numComp,3)] = f_100[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_101[i] + f_110[i]);                  } /* end of k0 loop */
                 } /* end of component loop i */  
1449              } /* end of k1 loop */              } /* end of k1 loop */
1450          } /* end of k2 loop */          } /* end of face 5 */
1451      } /* end of face 1 */          /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */
1452      if (m_faceOffset[2] > -1) {      } else {
1453          const double tmp0_2 = 0.044658198738520451079;          /* GENERATOR SNIP_INTERPOLATE_FACES TOP */
1454          const double tmp0_1 = 0.16666666666666666667;          if (m_faceOffset[0] > -1) {
1455          const double tmp0_0 = 0.62200846792814621559;              const double tmp0_2 = 0.044658198738520451079;
1456                const double tmp0_1 = 0.16666666666666666667;
1457                const double tmp0_0 = 0.62200846792814621559;
1458  #pragma omp parallel for  #pragma omp parallel for
1459          for (index_t k2=0; k2 < m_NE2; ++k2) {              for (index_t k2=0; k2 < m_NE2; ++k2) {
1460              for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1461                  const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));                      const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
1462                  const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));                      const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
1463                  const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));                      const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
1464                  const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));                      const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
1465                  double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                      double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1466                  for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1467                      o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_101[i]*tmp0_2 + tmp0_1*(f_001[i] + f_100[i]);                          o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_2 + tmp0_1*(f_001[i] + f_010[i]);
1468                      o[INDEX2(i,numComp,1)] = f_001[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_101[i]);                          o[INDEX2(i,numComp,1)] = f_001[i]*tmp0_2 + f_010[i]*tmp0_0 + tmp0_1*(f_000[i] + f_011[i]);
1469                      o[INDEX2(i,numComp,2)] = f_001[i]*tmp0_0 + f_100[i]*tmp0_2 + tmp0_1*(f_000[i] + f_101[i]);                          o[INDEX2(i,numComp,2)] = f_001[i]*tmp0_0 + f_010[i]*tmp0_2 + tmp0_1*(f_000[i] + f_011[i]);
1470                      o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_100[i]);                          o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_011[i]*tmp0_0 + tmp0_1*(f_001[i] + f_010[i]);
1471                  } /* end of component loop i */                      } /* end of component loop i */
1472              } /* end of k0 loop */                  } /* end of k1 loop */
1473          } /* end of k2 loop */              } /* end of k2 loop */
1474      } /* end of face 2 */          } /* end of face 0 */
1475      if (m_faceOffset[3] > -1) {          if (m_faceOffset[1] > -1) {
1476          const double tmp0_2 = 0.044658198738520451079;              const double tmp0_2 = 0.044658198738520451079;
1477          const double tmp0_1 = 0.62200846792814621559;              const double tmp0_1 = 0.62200846792814621559;
1478          const double tmp0_0 = 0.16666666666666666667;              const double tmp0_0 = 0.16666666666666666667;
1479  #pragma omp parallel for  #pragma omp parallel for
1480          for (index_t k2=0; k2 < m_NE2; ++k2) {              for (index_t k2=0; k2 < m_NE2; ++k2) {
1481              for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1482                  const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));                      const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
1483                  const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));                      const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
1484                  const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));                      const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
1485                  const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));                      const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
1486                  double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                      double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1487                  for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1488                      o[INDEX2(i,numComp,0)] = f_010[i]*tmp0_1 + f_111[i]*tmp0_2 + tmp0_0*(f_011[i] + f_110[i]);                          o[INDEX2(i,numComp,0)] = f_100[i]*tmp0_1 + f_111[i]*tmp0_2 + tmp0_0*(f_101[i] + f_110[i]);
1489                      o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_2 + f_110[i]*tmp0_1 + tmp0_0*(f_010[i] + f_111[i]);                          o[INDEX2(i,numComp,1)] = f_101[i]*tmp0_2 + f_110[i]*tmp0_1 + tmp0_0*(f_100[i] + f_111[i]);
1490                      o[INDEX2(i,numComp,2)] = f_011[i]*tmp0_1 + f_110[i]*tmp0_2 + tmp0_0*(f_010[i] + f_111[i]);                          o[INDEX2(i,numComp,2)] = f_101[i]*tmp0_1 + f_110[i]*tmp0_2 + tmp0_0*(f_100[i] + f_111[i]);
1491                      o[INDEX2(i,numComp,3)] = f_010[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_011[i] + f_110[i]);                          o[INDEX2(i,numComp,3)] = f_100[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_101[i] + f_110[i]);
1492                  } /* end of component loop i */                      } /* end of component loop i */
1493              } /* end of k0 loop */                  } /* end of k1 loop */
1494          } /* end of k2 loop */              } /* end of k2 loop */
1495      } /* end of face 3 */          } /* end of face 1 */
1496      if (m_faceOffset[4] > -1) {          if (m_faceOffset[2] > -1) {
1497          const double tmp0_2 = 0.044658198738520451079;              const double tmp0_2 = 0.044658198738520451079;
1498          const double tmp0_1 = 0.16666666666666666667;              const double tmp0_1 = 0.16666666666666666667;
1499          const double tmp0_0 = 0.62200846792814621559;              const double tmp0_0 = 0.62200846792814621559;
1500  #pragma omp parallel for  #pragma omp parallel for
1501          for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k2=0; k2 < m_NE2; ++k2) {
1502              for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1503                  const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));                      const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
1504                  const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));                      const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
1505                  const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));                      const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
1506                  const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));                      const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
1507                  double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                      double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1508                  for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1509                      o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_110[i]*tmp0_2 + tmp0_1*(f_010[i] + f_100[i]);                          o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_101[i]*tmp0_2 + tmp0_1*(f_001[i] + f_100[i]);
1510                      o[INDEX2(i,numComp,1)] = f_010[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_110[i]);                          o[INDEX2(i,numComp,1)] = f_001[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_101[i]);
1511                      o[INDEX2(i,numComp,2)] = f_010[i]*tmp0_0 + f_100[i]*tmp0_2 + tmp0_1*(f_000[i] + f_110[i]);                          o[INDEX2(i,numComp,2)] = f_001[i]*tmp0_0 + f_100[i]*tmp0_2 + tmp0_1*(f_000[i] + f_101[i]);
1512                      o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_110[i]*tmp0_0 + tmp0_1*(f_010[i] + f_100[i]);                          o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_100[i]);
1513                  } /* end of component loop i */                      } /* end of component loop i */
1514              } /* end of k0 loop */                  } /* end of k0 loop */
1515          } /* end of k1 loop */              } /* end of k2 loop */
1516      } /* end of face 4 */          } /* end of face 2 */
1517      if (m_faceOffset[5] > -1) {          if (m_faceOffset[3] > -1) {
1518          const double tmp0_2 = 0.044658198738520451079;              const double tmp0_2 = 0.044658198738520451079;
1519          const double tmp0_1 = 0.16666666666666666667;              const double tmp0_1 = 0.62200846792814621559;
1520          const double tmp0_0 = 0.62200846792814621559;              const double tmp0_0 = 0.16666666666666666667;
1521  #pragma omp parallel for  #pragma omp parallel for
1522          for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k2=0; k2 < m_NE2; ++k2) {
1523              for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1524                  const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));                      const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
1525                  const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));                      const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
1526                  const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));                      const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
1527                  const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));                      const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
1528                  double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                      double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1529                  for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1530                      o[INDEX2(i,numComp,0)] = f_001[i]*tmp0_0 + f_111[i]*tmp0_2 + tmp0_1*(f_011[i] + f_101[i]);                          o[INDEX2(i,numComp,0)] = f_010[i]*tmp0_1 + f_111[i]*tmp0_2 + tmp0_0*(f_011[i] + f_110[i]);
1531                      o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_111[i]);                          o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_2 + f_110[i]*tmp0_1 + tmp0_0*(f_010[i] + f_111[i]);
1532                      o[INDEX2(i,numComp,2)] = f_011[i]*tmp0_0 + f_101[i]*tmp0_2 + tmp0_1*(f_001[i] + f_111[i]);                          o[INDEX2(i,numComp,2)] = f_011[i]*tmp0_1 + f_110[i]*tmp0_2 + tmp0_0*(f_010[i] + f_111[i]);
1533                      o[INDEX2(i,numComp,3)] = f_001[i]*tmp0_2 + f_111[i]*tmp0_0 + tmp0_1*(f_011[i] + f_101[i]);                          o[INDEX2(i,numComp,3)] = f_010[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_011[i] + f_110[i]);
1534                  } /* end of component loop i */                      } /* end of component loop i */
1535              } /* end of k0 loop */                  } /* end of k0 loop */
1536          } /* end of k1 loop */              } /* end of k2 loop */
1537      } /* end of face 5 */          } /* end of face 3 */
1538      /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */          if (m_faceOffset[4] > -1) {
1539                const double tmp0_2 = 0.044658198738520451079;
1540                const double tmp0_1 = 0.16666666666666666667;
1541                const double tmp0_0 = 0.62200846792814621559;
1542    #pragma omp parallel for
1543                for (index_t k1=0; k1 < m_NE1; ++k1) {
1544                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1545                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
1546                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
1547                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
1548                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
1549                        double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1550                        for (index_t i=0; i < numComp; ++i) {
1551                            o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_110[i]*tmp0_2 + tmp0_1*(f_010[i] + f_100[i]);
1552                            o[INDEX2(i,numComp,1)] = f_010[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_110[i]);
1553                            o[INDEX2(i,numComp,2)] = f_010[i]*tmp0_0 + f_100[i]*tmp0_2 + tmp0_1*(f_000[i] + f_110[i]);
1554                            o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_110[i]*tmp0_0 + tmp0_1*(f_010[i] + f_100[i]);
1555                        } /* end of component loop i */
1556                    } /* end of k0 loop */
1557                } /* end of k1 loop */
1558            } /* end of face 4 */
1559            if (m_faceOffset[5] > -1) {
1560                const double tmp0_2 = 0.044658198738520451079;
1561                const double tmp0_1 = 0.16666666666666666667;
1562                const double tmp0_0 = 0.62200846792814621559;
1563    #pragma omp parallel for
1564                for (index_t k1=0; k1 < m_NE1; ++k1) {
1565                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1566                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
1567                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
1568                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
1569                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
1570                        double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1571                        for (index_t i=0; i < numComp; ++i) {
1572                            o[INDEX2(i,numComp,0)] = f_001[i]*tmp0_0 + f_111[i]*tmp0_2 + tmp0_1*(f_011[i] + f_101[i]);
1573                            o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_111[i]);
1574                            o[INDEX2(i,numComp,2)] = f_011[i]*tmp0_0 + f_101[i]*tmp0_2 + tmp0_1*(f_001[i] + f_111[i]);
1575                            o[INDEX2(i,numComp,3)] = f_001[i]*tmp0_2 + f_111[i]*tmp0_0 + tmp0_1*(f_011[i] + f_101[i]);
1576                        } /* end of component loop i */
1577                    } /* end of k0 loop */
1578                } /* end of k1 loop */
1579            } /* end of face 5 */
1580            /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
1581        }
1582  }  }
1583    
1584  } // end of namespace ripley  } // end of namespace ripley

Legend:
Removed from v.3710  
changed lines
  Added in v.3711

  ViewVC Help
Powered by ViewVC 1.1.26