/[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 3722 by caltinay, Wed Dec 7 05:53:22 2011 UTC revision 3731 by caltinay, Fri Dec 9 01:22:44 2011 UTC
# Line 268  void Brick::setToGradient(escript::Data& Line 268  void Brick::setToGradient(escript::Data&
268      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
269      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
270      const double h2 = m_l1/m_gNE2;      const double h2 = m_l1/m_gNE2;
271        const double C0 = .044658198738520451079;
272        const double C1 = .16666666666666666667;
273        const double C2 = .21132486540518711775;
274        const double C3 = .25;
275        const double C4 = .5;
276        const double C5 = .62200846792814621559;
277        const double C6 = .78867513459481288225;
278    
279      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
280          /* GENERATOR SNIP_GRAD_ELEMENTS TOP */          /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */
         const double tmp0_22 = -0.044658198738520451079/h1;  
         const double tmp0_16 = 0.16666666666666666667/h0;  
         const double tmp0_33 = -0.62200846792814621559/h1;  
         const double tmp0_0 = -0.62200846792814621559/h0;  
         const double tmp0_21 = -0.16666666666666666667/h1;  
         const double tmp0_17 = 0.62200846792814621559/h0;  
         const double tmp0_52 = -0.044658198738520451079/h2;  
         const double tmp0_1 = -0.16666666666666666667/h0;  
         const double tmp0_20 = -0.62200846792814621559/h1;  
         const double tmp0_14 = -0.044658198738520451079/h0;  
         const double tmp0_53 = -0.62200846792814621559/h2;  
         const double tmp0_49 = 0.16666666666666666667/h2;  
         const double tmp0_2 = 0.16666666666666666667/h0;  
         const double tmp0_27 = -0.044658198738520451079/h1;  
         const double tmp0_15 = -0.16666666666666666667/h0;  
         const double tmp0_50 = -0.16666666666666666667/h2;  
         const double tmp0_48 = 0.62200846792814621559/h2;  
         const double tmp0_3 = 0.044658198738520451079/h0;  
         const double tmp0_26 = -0.16666666666666666667/h1;  
         const double tmp0_12 = -0.62200846792814621559/h0;  
         const double tmp0_51 = 0.044658198738520451079/h2;  
         const double tmp0_25 = 0.62200846792814621559/h1;  
         const double tmp0_13 = 0.16666666666666666667/h0;  
         const double tmp0_56 = 0.16666666666666666667/h2;  
         const double tmp0_24 = 0.16666666666666666667/h1;  
         const double tmp0_10 = 0.62200846792814621559/h0;  
         const double tmp0_57 = 0.62200846792814621559/h2;  
         const double tmp0_11 = -0.16666666666666666667/h0;  
         const double tmp0_54 = -0.044658198738520451079/h2;  
         const double tmp0_38 = 0.16666666666666666667/h1;  
         const double tmp0_34 = -0.044658198738520451079/h1;  
         const double tmp0_42 = 0.16666666666666666667/h2;  
         const double tmp0_35 = -0.16666666666666666667/h1;  
         const double tmp0_36 = -0.62200846792814621559/h1;  
         const double tmp0_41 = 0.62200846792814621559/h2;  
         const double tmp0_8 = 0.044658198738520451079/h0;  
         const double tmp0_37 = 0.62200846792814621559/h1;  
         const double tmp0_29 = 0.16666666666666666667/h1;  
         const double tmp0_40 = -0.62200846792814621559/h2;  
         const double tmp0_9 = 0.16666666666666666667/h0;  
         const double tmp0_30 = 0.62200846792814621559/h1;  
         const double tmp0_28 = -0.16666666666666666667/h1;  
         const double tmp0_43 = 0.044658198738520451079/h2;  
         const double tmp0_32 = 0.16666666666666666667/h1;  
         const double tmp0_31 = 0.044658198738520451079/h1;  
         const double tmp0_39 = 0.044658198738520451079/h1;  
         const double tmp0_58 = -0.62200846792814621559/h2;  
         const double tmp0_55 = 0.044658198738520451079/h2;  
         const double tmp0_18 = -0.62200846792814621559/h0;  
         const double tmp0_45 = -0.16666666666666666667/h2;  
         const double tmp0_59 = -0.16666666666666666667/h2;  
         const double tmp0_4 = -0.044658198738520451079/h0;  
         const double tmp0_19 = 0.044658198738520451079/h0;  
         const double tmp0_44 = -0.044658198738520451079/h2;  
         const double tmp0_5 = 0.62200846792814621559/h0;  
         const double tmp0_47 = 0.16666666666666666667/h2;  
         const double tmp0_6 = -0.16666666666666666667/h0;  
         const double tmp0_23 = 0.044658198738520451079/h1;  
         const double tmp0_46 = -0.16666666666666666667/h2;  
         const double tmp0_7 = -0.044658198738520451079/h0;  
281  #pragma omp parallel for  #pragma omp parallel for
282          for (index_t k2=0; k2 < m_NE2; ++k2) {          for (index_t k2=0; k2 < m_NE2; ++k2) {
283              for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 344  void Brick::setToGradient(escript::Data& Line 292  void Brick::setToGradient(escript::Data&
292                      const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));                      const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
293                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
294                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
295                          o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);                          const double V0=((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
296                          o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_20 + f_010[i]*tmp0_25 + f_101[i]*tmp0_22 + f_111[i]*tmp0_23 + tmp0_21*(f_001[i] + f_100[i]) + tmp0_24*(f_011[i] + f_110[i]);                          const double V1=((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
297                          o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_40 + f_001[i]*tmp0_41 + f_110[i]*tmp0_44 + f_111[i]*tmp0_43 + tmp0_42*(f_011[i] + f_101[i]) + tmp0_45*(f_010[i] + f_100[i]);                          const double V2=((f_101[i]-f_001[i])*C5 + (f_110[i]-f_010[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
298                          o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);                          const double V3=((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
299                          o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_27 + f_010[i]*tmp0_32 + f_011[i]*tmp0_31 + f_100[i]*tmp0_33 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;                          const double V4=((f_010[i]-f_000[i])*C5 + (f_111[i]-f_101[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
300                          o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_52 + f_011[i]*tmp0_51 + f_100[i]*tmp0_53 + f_101[i]*tmp0_48 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;                          const double V5=((f_110[i]-f_100[i])*C5 + (f_011[i]-f_001[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
301                          o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;                          const double V6=((f_011[i]-f_001[i])*C5 + (f_110[i]-f_100[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
302                          o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_20 + f_010[i]*tmp0_25 + f_101[i]*tmp0_22 + f_111[i]*tmp0_23 + tmp0_21*(f_001[i] + f_100[i]) + tmp0_24*(f_011[i] + f_110[i]);                          const double V7=((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
303                          o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_53 + f_011[i]*tmp0_48 + f_100[i]*tmp0_52 + f_101[i]*tmp0_51 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;                          const double V8=((f_001[i]-f_000[i])*C5 + (f_111[i]-f_110[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
304                          o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;                          const double V9=((f_101[i]-f_100[i])*C5 + (f_011[i]-f_010[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
305                          o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_27 + f_010[i]*tmp0_32 + f_011[i]*tmp0_31 + f_100[i]*tmp0_33 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;                          const double V10=((f_011[i]-f_010[i])*C5 + (f_101[i]-f_100[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
306                          o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_54 + f_001[i]*tmp0_55 + f_110[i]*tmp0_58 + f_111[i]*tmp0_57 + tmp0_56*(f_011[i] + f_101[i]) + tmp0_59*(f_010[i] + f_100[i]);                          const double V11=((f_111[i]-f_110[i])*C5 + (f_001[i]-f_000[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
307                          o[INDEX3(i,0,4,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;                          o[INDEX3(i,0,0,numComp,3)] = V0;
308                          o[INDEX3(i,1,4,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_33 + f_010[i]*tmp0_32 + f_011[i]*tmp0_30 + f_100[i]*tmp0_27 + f_101[i]*tmp0_28 + f_110[i]*tmp0_31 + f_111[i]*tmp0_29;                          o[INDEX3(i,1,0,numComp,3)] = V4;
309                          o[INDEX3(i,2,4,numComp,3)] = f_000[i]*tmp0_40 + f_001[i]*tmp0_41 + f_110[i]*tmp0_44 + f_111[i]*tmp0_43 + tmp0_42*(f_011[i] + f_101[i]) + tmp0_45*(f_010[i] + f_100[i]);                          o[INDEX3(i,2,0,numComp,3)] = V8;
310                          o[INDEX3(i,0,5,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;                          o[INDEX3(i,0,1,numComp,3)] = V0;
311                          o[INDEX3(i,1,5,numComp,3)] = f_000[i]*tmp0_34 + f_010[i]*tmp0_39 + f_101[i]*tmp0_36 + f_111[i]*tmp0_37 + tmp0_35*(f_001[i] + f_100[i]) + tmp0_38*(f_011[i] + f_110[i]);                          o[INDEX3(i,1,1,numComp,3)] = V5;
312                          o[INDEX3(i,2,5,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_52 + f_011[i]*tmp0_51 + f_100[i]*tmp0_53 + f_101[i]*tmp0_48 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;                          o[INDEX3(i,2,1,numComp,3)] = V9;
313                          o[INDEX3(i,0,6,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);                          o[INDEX3(i,0,2,numComp,3)] = V1;
314                          o[INDEX3(i,1,6,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_33 + f_010[i]*tmp0_32 + f_011[i]*tmp0_30 + f_100[i]*tmp0_27 + f_101[i]*tmp0_28 + f_110[i]*tmp0_31 + f_111[i]*tmp0_29;                          o[INDEX3(i,1,2,numComp,3)] = V4;
315                          o[INDEX3(i,2,6,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_53 + f_011[i]*tmp0_48 + f_100[i]*tmp0_52 + f_101[i]*tmp0_51 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;                          o[INDEX3(i,2,2,numComp,3)] = V10;
316                          o[INDEX3(i,0,7,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);                          o[INDEX3(i,0,3,numComp,3)] = V1;
317                          o[INDEX3(i,1,7,numComp,3)] = f_000[i]*tmp0_34 + f_010[i]*tmp0_39 + f_101[i]*tmp0_36 + f_111[i]*tmp0_37 + tmp0_35*(f_001[i] + f_100[i]) + tmp0_38*(f_011[i] + f_110[i]);                          o[INDEX3(i,1,3,numComp,3)] = V5;
318                          o[INDEX3(i,2,7,numComp,3)] = f_000[i]*tmp0_54 + f_001[i]*tmp0_55 + f_110[i]*tmp0_58 + f_111[i]*tmp0_57 + tmp0_56*(f_011[i] + f_101[i]) + tmp0_59*(f_010[i] + f_100[i]);                          o[INDEX3(i,2,3,numComp,3)] = V11;
319                            o[INDEX3(i,0,4,numComp,3)] = V2;
320                            o[INDEX3(i,1,4,numComp,3)] = V6;
321                            o[INDEX3(i,2,4,numComp,3)] = V8;
322                            o[INDEX3(i,0,5,numComp,3)] = V2;
323                            o[INDEX3(i,1,5,numComp,3)] = V7;
324                            o[INDEX3(i,2,5,numComp,3)] = V9;
325                            o[INDEX3(i,0,6,numComp,3)] = V3;
326                            o[INDEX3(i,1,6,numComp,3)] = V6;
327                            o[INDEX3(i,2,6,numComp,3)] = V10;
328                            o[INDEX3(i,0,7,numComp,3)] = V3;
329                            o[INDEX3(i,1,7,numComp,3)] = V7;
330                            o[INDEX3(i,2,7,numComp,3)] = V11;
331                      } /* end of component loop i */                      } /* end of component loop i */
332                  } /* end of k0 loop */                  } /* end of k0 loop */
333              } /* end of k1 loop */              } /* end of k1 loop */
334          } /* end of k2 loop */          } /* end of k2 loop */
335          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
336      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
337          /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */          /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */
         const double tmp0_0 = -0.25/h0;  
         const double tmp0_1 = 0.25/h0;  
         const double tmp0_2 = -0.25/h1;  
         const double tmp0_3 = 0.25/h1;  
         const double tmp0_4 = -0.25/h2;  
         const double tmp0_5 = 0.25/h2;  
338  #pragma omp parallel for  #pragma omp parallel for
339          for (index_t k2=0; k2 < m_NE2; ++k2) {          for (index_t k2=0; k2 < m_NE2; ++k2) {
340              for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 395  void Brick::setToGradient(escript::Data& Line 349  void Brick::setToGradient(escript::Data&
349                      const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));                      const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
350                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
351                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
352                          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]);                          o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;
353                          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]);                          o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;
354                          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]);                          o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / h2;
355                      } /* end of component loop i */                      } /* end of component loop i */
356                  } /* end of k0 loop */                  } /* end of k0 loop */
357              } /* end of k1 loop */              } /* end of k1 loop */
358          } /* end of k2 loop */          } /* end of k2 loop */
359          /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */          /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
360      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
361          /* GENERATOR SNIP_GRAD_FACES TOP */          /*** GENERATOR SNIP_GRAD_FACES TOP */
362  #pragma omp parallel  #pragma omp parallel
363          {          {
364              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
                 const double tmp0_22 = 0.21132486540518711775/h1;  
                 const double tmp0_16 = 0.16666666666666666667/h0;  
                 const double tmp0_33 = 0.21132486540518711775/h2;  
                 const double tmp0_0 = -0.62200846792814621559/h0;  
                 const double tmp0_21 = -0.21132486540518711775/h1;  
                 const double tmp0_17 = 0.62200846792814621559/h0;  
                 const double tmp0_1 = -0.16666666666666666667/h0;  
                 const double tmp0_20 = -0.78867513459481288225/h1;  
                 const double tmp0_14 = -0.044658198738520451079/h0;  
                 const double tmp0_2 = 0.16666666666666666667/h0;  
                 const double tmp0_27 = 0.21132486540518711775/h1;  
                 const double tmp0_15 = -0.16666666666666666667/h0;  
                 const double tmp0_3 = 0.044658198738520451079/h0;  
                 const double tmp0_26 = 0.78867513459481288225/h1;  
                 const double tmp0_12 = -0.62200846792814621559/h0;  
                 const double tmp0_25 = -0.78867513459481288225/h1;  
                 const double tmp0_13 = 0.16666666666666666667/h0;  
                 const double tmp0_24 = -0.21132486540518711775/h1;  
                 const double tmp0_10 = 0.62200846792814621559/h0;  
                 const double tmp0_11 = -0.16666666666666666667/h0;  
                 const double tmp0_34 = 0.78867513459481288225/h2;  
                 const double tmp0_35 = -0.78867513459481288225/h2;  
                 const double tmp0_8 = 0.044658198738520451079/h0;  
                 const double tmp0_29 = 0.78867513459481288225/h2;  
                 const double tmp0_9 = 0.16666666666666666667/h0;  
                 const double tmp0_30 = 0.21132486540518711775/h2;  
                 const double tmp0_28 = -0.78867513459481288225/h2;  
                 const double tmp0_32 = -0.21132486540518711775/h2;  
                 const double tmp0_31 = -0.21132486540518711775/h2;  
                 const double tmp0_18 = -0.62200846792814621559/h0;  
                 const double tmp0_4 = -0.044658198738520451079/h0;  
                 const double tmp0_19 = 0.044658198738520451079/h0;  
                 const double tmp0_5 = 0.62200846792814621559/h0;  
                 const double tmp0_6 = -0.16666666666666666667/h0;  
                 const double tmp0_23 = 0.78867513459481288225/h1;  
                 const double tmp0_7 = -0.044658198738520451079/h0;  
365  #pragma omp for nowait  #pragma omp for nowait
366                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
367                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 457  void Brick::setToGradient(escript::Data& Line 375  void Brick::setToGradient(escript::Data&
375                          const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));                          const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
376                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
377                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
378                              o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);                              const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / h1;
379                              o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_20 + f_001[i]*tmp0_21 + f_010[i]*tmp0_23 + f_011[i]*tmp0_22;                              const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / h1;
380                              o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_28 + f_001[i]*tmp0_29 + f_010[i]*tmp0_31 + f_011[i]*tmp0_30;                              const double V2=((f_001[i]-f_000[i])*C6 + (f_010[i]-f_011[i])*C2) / h2;
381                              o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;                              const double V3=((f_001[i]-f_000[i])*C2 + (f_011[i]-f_010[i])*C6) / h2;
382                              o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_20 + f_001[i]*tmp0_21 + f_010[i]*tmp0_23 + f_011[i]*tmp0_22;                              o[INDEX3(i,0,0,numComp,3)] = ((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
383                              o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_32 + f_001[i]*tmp0_33 + f_010[i]*tmp0_35 + f_011[i]*tmp0_34;                              o[INDEX3(i,1,0,numComp,3)] = V0;
384                              o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;                              o[INDEX3(i,2,0,numComp,3)] = V2;
385                              o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_24 + f_001[i]*tmp0_25 + f_010[i]*tmp0_27 + f_011[i]*tmp0_26;                              o[INDEX3(i,0,1,numComp,3)] = ((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
386                              o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_28 + f_001[i]*tmp0_29 + f_010[i]*tmp0_31 + f_011[i]*tmp0_30;                              o[INDEX3(i,1,1,numComp,3)] = V0;
387                              o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);                              o[INDEX3(i,2,1,numComp,3)] = V3;
388                              o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_24 + f_001[i]*tmp0_25 + f_010[i]*tmp0_27 + f_011[i]*tmp0_26;                              o[INDEX3(i,0,2,numComp,3)] = ((f_101[i]-f_001[i])*C5 + (f_110[i]-f_010[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
389                              o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_32 + f_001[i]*tmp0_33 + f_010[i]*tmp0_35 + f_011[i]*tmp0_34;                              o[INDEX3(i,1,2,numComp,3)] = V1;
390                                o[INDEX3(i,2,2,numComp,3)] = V2;
391                                o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
392                                o[INDEX3(i,1,3,numComp,3)] = V1;
393                                o[INDEX3(i,2,3,numComp,3)] = V3;
394                          } /* end of component loop i */                          } /* end of component loop i */
395                      } /* end of k1 loop */                      } /* end of k1 loop */
396                  } /* end of k2 loop */                  } /* end of k2 loop */
397              } /* end of face 0 */              } /* end of face 0 */
398              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
                 const double tmp0_22 = 0.78867513459481288225/h1;  
                 const double tmp0_16 = 0.16666666666666666667/h0;  
                 const double tmp0_33 = 0.78867513459481288225/h2;  
                 const double tmp0_0 = -0.62200846792814621559/h0;  
                 const double tmp0_21 = 0.21132486540518711775/h1;  
                 const double tmp0_17 = 0.62200846792814621559/h0;  
                 const double tmp0_1 = -0.16666666666666666667/h0;  
                 const double tmp0_20 = -0.21132486540518711775/h1;  
                 const double tmp0_14 = -0.044658198738520451079/h0;  
                 const double tmp0_2 = 0.16666666666666666667/h0;  
                 const double tmp0_27 = -0.21132486540518711775/h1;  
                 const double tmp0_15 = -0.16666666666666666667/h0;  
                 const double tmp0_3 = 0.044658198738520451079/h0;  
                 const double tmp0_26 = 0.21132486540518711775/h1;  
                 const double tmp0_12 = -0.62200846792814621559/h0;  
                 const double tmp0_25 = 0.78867513459481288225/h1;  
                 const double tmp0_13 = 0.16666666666666666667/h0;  
                 const double tmp0_24 = -0.78867513459481288225/h1;  
                 const double tmp0_10 = 0.62200846792814621559/h0;  
                 const double tmp0_11 = -0.16666666666666666667/h0;  
                 const double tmp0_34 = -0.78867513459481288225/h2;  
                 const double tmp0_35 = -0.21132486540518711775/h2;  
                 const double tmp0_8 = 0.044658198738520451079/h0;  
                 const double tmp0_29 = 0.21132486540518711775/h2;  
                 const double tmp0_9 = 0.16666666666666666667/h0;  
                 const double tmp0_30 = -0.21132486540518711775/h2;  
                 const double tmp0_28 = 0.78867513459481288225/h2;  
                 const double tmp0_32 = 0.21132486540518711775/h2;  
                 const double tmp0_31 = -0.78867513459481288225/h2;  
                 const double tmp0_18 = -0.62200846792814621559/h0;  
                 const double tmp0_4 = -0.044658198738520451079/h0;  
                 const double tmp0_19 = 0.044658198738520451079/h0;  
                 const double tmp0_5 = 0.62200846792814621559/h0;  
                 const double tmp0_6 = -0.16666666666666666667/h0;  
                 const double tmp0_23 = -0.78867513459481288225/h1;  
                 const double tmp0_7 = -0.044658198738520451079/h0;  
399  #pragma omp for nowait  #pragma omp for nowait
400                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
401                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 523  void Brick::setToGradient(escript::Data& Line 409  void Brick::setToGradient(escript::Data&
409                          const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));                          const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
410                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
411                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
412                              o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);                              const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;
413                              o[INDEX3(i,1,0,numComp,3)] = f_100[i]*tmp0_23 + f_101[i]*tmp0_20 + f_110[i]*tmp0_22 + f_111[i]*tmp0_21;                              const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;
414                              o[INDEX3(i,2,0,numComp,3)] = f_100[i]*tmp0_31 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;                              const double V2=((f_101[i]-f_100[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;
415                              o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;                              const double V3=((f_101[i]-f_100[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;
416                              o[INDEX3(i,1,1,numComp,3)] = f_100[i]*tmp0_23 + f_101[i]*tmp0_20 + f_110[i]*tmp0_22 + f_111[i]*tmp0_21;                              o[INDEX3(i,0,0,numComp,3)] = ((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
417                              o[INDEX3(i,2,1,numComp,3)] = f_100[i]*tmp0_35 + f_101[i]*tmp0_32 + f_110[i]*tmp0_34 + f_111[i]*tmp0_33;                              o[INDEX3(i,1,0,numComp,3)] = V0;
418                              o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;                              o[INDEX3(i,2,0,numComp,3)] = V2;
419                              o[INDEX3(i,1,2,numComp,3)] = f_100[i]*tmp0_27 + f_101[i]*tmp0_24 + f_110[i]*tmp0_26 + f_111[i]*tmp0_25;                              o[INDEX3(i,0,1,numComp,3)] = ((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
420                              o[INDEX3(i,2,2,numComp,3)] = f_100[i]*tmp0_31 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;                              o[INDEX3(i,1,1,numComp,3)] = V0;
421                              o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);                              o[INDEX3(i,2,1,numComp,3)] = V3;
422                              o[INDEX3(i,1,3,numComp,3)] = f_100[i]*tmp0_27 + f_101[i]*tmp0_24 + f_110[i]*tmp0_26 + f_111[i]*tmp0_25;                              o[INDEX3(i,0,2,numComp,3)] = ((f_101[i]-f_001[i])*C5 + (f_110[i]-f_010[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
423                              o[INDEX3(i,2,3,numComp,3)] = f_100[i]*tmp0_35 + f_101[i]*tmp0_32 + f_110[i]*tmp0_34 + f_111[i]*tmp0_33;                              o[INDEX3(i,1,2,numComp,3)] = V1;
424                                o[INDEX3(i,2,2,numComp,3)] = V2;
425                                o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
426                                o[INDEX3(i,1,3,numComp,3)] = V1;
427                                o[INDEX3(i,2,3,numComp,3)] = V3;
428                          } /* end of component loop i */                          } /* end of component loop i */
429                      } /* end of k1 loop */                      } /* end of k1 loop */
430                  } /* end of k2 loop */                  } /* end of k2 loop */
431              } /* end of face 1 */              } /* end of face 1 */
432              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
                 const double tmp0_22 = -0.044658198738520451079/h1;  
                 const double tmp0_16 = -0.16666666666666666667/h1;  
                 const double tmp0_33 = 0.21132486540518711775/h2;  
                 const double tmp0_0 = -0.78867513459481288225/h0;  
                 const double tmp0_21 = 0.16666666666666666667/h1;  
                 const double tmp0_17 = -0.62200846792814621559/h1;  
                 const double tmp0_1 = -0.21132486540518711775/h0;  
                 const double tmp0_20 = 0.044658198738520451079/h1;  
                 const double tmp0_14 = -0.16666666666666666667/h1;  
                 const double tmp0_2 = 0.21132486540518711775/h0;  
                 const double tmp0_27 = 0.044658198738520451079/h1;  
                 const double tmp0_15 = -0.044658198738520451079/h1;  
                 const double tmp0_3 = 0.78867513459481288225/h0;  
                 const double tmp0_26 = 0.16666666666666666667/h1;  
                 const double tmp0_12 = 0.16666666666666666667/h1;  
                 const double tmp0_25 = 0.62200846792814621559/h1;  
                 const double tmp0_13 = 0.62200846792814621559/h1;  
                 const double tmp0_24 = -0.62200846792814621559/h1;  
                 const double tmp0_10 = -0.044658198738520451079/h1;  
                 const double tmp0_11 = 0.044658198738520451079/h1;  
                 const double tmp0_34 = 0.78867513459481288225/h2;  
                 const double tmp0_35 = -0.78867513459481288225/h2;  
                 const double tmp0_8 = -0.62200846792814621559/h1;  
                 const double tmp0_29 = 0.78867513459481288225/h2;  
                 const double tmp0_9 = -0.16666666666666666667/h1;  
                 const double tmp0_30 = 0.21132486540518711775/h2;  
                 const double tmp0_28 = -0.78867513459481288225/h2;  
                 const double tmp0_32 = -0.21132486540518711775/h2;  
                 const double tmp0_31 = -0.21132486540518711775/h2;  
                 const double tmp0_18 = 0.16666666666666666667/h1;  
                 const double tmp0_4 = -0.21132486540518711775/h0;  
                 const double tmp0_19 = 0.62200846792814621559/h1;  
                 const double tmp0_5 = -0.78867513459481288225/h0;  
                 const double tmp0_6 = 0.78867513459481288225/h0;  
                 const double tmp0_23 = -0.16666666666666666667/h1;  
                 const double tmp0_7 = 0.21132486540518711775/h0;  
433  #pragma omp for nowait  #pragma omp for nowait
434                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
435                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 589  void Brick::setToGradient(escript::Data& Line 443  void Brick::setToGradient(escript::Data&
443                          const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));                          const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
444                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
445                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
446                              o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_001[i]*tmp0_1 + f_100[i]*tmp0_3 + f_101[i]*tmp0_2;                              const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / h0;
447                              o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_8 + f_010[i]*tmp0_13 + f_101[i]*tmp0_10 + f_111[i]*tmp0_11 + tmp0_12*(f_011[i] + f_110[i]) + tmp0_9*(f_001[i] + f_100[i]);                              const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / h2;
448                              o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_28 + f_001[i]*tmp0_29 + f_100[i]*tmp0_31 + f_101[i]*tmp0_30;                              const double V2=((f_001[i]-f_000[i])*C2 + (f_101[i]-f_100[i])*C6) / h2;
449                              o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_0 + f_001[i]*tmp0_1 + f_100[i]*tmp0_3 + f_101[i]*tmp0_2;                              o[INDEX3(i,0,0,numComp,3)] = V0;
450                              o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_14 + f_001[i]*tmp0_15 + f_010[i]*tmp0_21 + f_011[i]*tmp0_20 + f_100[i]*tmp0_17 + f_101[i]*tmp0_16 + f_110[i]*tmp0_19 + f_111[i]*tmp0_18;                              o[INDEX3(i,1,0,numComp,3)] = ((f_010[i]-f_000[i])*C5 + (f_111[i]-f_101[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
451                              o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_32 + f_001[i]*tmp0_33 + f_100[i]*tmp0_35 + f_101[i]*tmp0_34;                              o[INDEX3(i,2,0,numComp,3)] = V1;
452                              o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_4 + f_001[i]*tmp0_5 + f_100[i]*tmp0_7 + f_101[i]*tmp0_6;                              o[INDEX3(i,0,1,numComp,3)] = V0;
453                              o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_14 + f_001[i]*tmp0_17 + f_010[i]*tmp0_21 + f_011[i]*tmp0_19 + f_100[i]*tmp0_15 + f_101[i]*tmp0_16 + f_110[i]*tmp0_20 + f_111[i]*tmp0_18;                              o[INDEX3(i,1,1,numComp,3)] = ((f_110[i]-f_100[i])*C5 + (f_011[i]-f_001[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
454                              o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_28 + f_001[i]*tmp0_29 + f_100[i]*tmp0_31 + f_101[i]*tmp0_30;                              o[INDEX3(i,2,1,numComp,3)] = V2;
455                              o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_4 + f_001[i]*tmp0_5 + f_100[i]*tmp0_7 + f_101[i]*tmp0_6;                              o[INDEX3(i,0,2,numComp,3)] = V0;
456                              o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_22 + f_010[i]*tmp0_27 + f_101[i]*tmp0_24 + f_111[i]*tmp0_25 + tmp0_23*(f_001[i] + f_100[i]) + tmp0_26*(f_011[i] + f_110[i]);                              o[INDEX3(i,1,2,numComp,3)] = ((f_011[i]-f_001[i])*C5 + (f_110[i]-f_100[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
457                              o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_32 + f_001[i]*tmp0_33 + f_100[i]*tmp0_35 + f_101[i]*tmp0_34;                              o[INDEX3(i,2,2,numComp,3)] = V1;
458                                o[INDEX3(i,0,3,numComp,3)] = V0;
459                                o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
460                                o[INDEX3(i,2,3,numComp,3)] = V2;
461                          } /* end of component loop i */                          } /* end of component loop i */
462                      } /* end of k0 loop */                      } /* end of k0 loop */
463                  } /* end of k2 loop */                  } /* end of k2 loop */
464              } /* end of face 2 */              } /* end of face 2 */
465              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
                 const double tmp0_22 = 0.16666666666666666667/h1;  
                 const double tmp0_16 = 0.16666666666666666667/h1;  
                 const double tmp0_33 = -0.78867513459481288225/h2;  
                 const double tmp0_0 = -0.21132486540518711775/h0;  
                 const double tmp0_21 = -0.62200846792814621559/h1;  
                 const double tmp0_17 = 0.16666666666666666667/h1;  
                 const double tmp0_1 = 0.78867513459481288225/h0;  
                 const double tmp0_20 = -0.16666666666666666667/h1;  
                 const double tmp0_14 = 0.044658198738520451079/h1;  
                 const double tmp0_2 = -0.78867513459481288225/h0;  
                 const double tmp0_27 = -0.62200846792814621559/h1;  
                 const double tmp0_15 = 0.62200846792814621559/h1;  
                 const double tmp0_3 = 0.21132486540518711775/h0;  
                 const double tmp0_26 = -0.16666666666666666667/h1;  
                 const double tmp0_12 = -0.16666666666666666667/h1;  
                 const double tmp0_25 = -0.044658198738520451079/h1;  
                 const double tmp0_13 = -0.044658198738520451079/h1;  
                 const double tmp0_24 = 0.62200846792814621559/h1;  
                 const double tmp0_10 = 0.044658198738520451079/h1;  
                 const double tmp0_11 = -0.62200846792814621559/h1;  
                 const double tmp0_34 = -0.21132486540518711775/h2;  
                 const double tmp0_35 = 0.78867513459481288225/h2;  
                 const double tmp0_8 = 0.16666666666666666667/h1;  
                 const double tmp0_29 = -0.21132486540518711775/h2;  
                 const double tmp0_9 = 0.62200846792814621559/h1;  
                 const double tmp0_30 = -0.78867513459481288225/h2;  
                 const double tmp0_28 = 0.78867513459481288225/h2;  
                 const double tmp0_32 = 0.21132486540518711775/h2;  
                 const double tmp0_31 = 0.21132486540518711775/h2;  
                 const double tmp0_18 = -0.16666666666666666667/h1;  
                 const double tmp0_4 = -0.78867513459481288225/h0;  
                 const double tmp0_19 = -0.044658198738520451079/h1;  
                 const double tmp0_5 = 0.21132486540518711775/h0;  
                 const double tmp0_6 = -0.21132486540518711775/h0;  
                 const double tmp0_23 = 0.044658198738520451079/h1;  
                 const double tmp0_7 = 0.78867513459481288225/h0;  
466  #pragma omp for nowait  #pragma omp for nowait
467                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
468                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 655  void Brick::setToGradient(escript::Data& Line 476  void Brick::setToGradient(escript::Data&
476                          const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));                          const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
477                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
478                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
479                              o[INDEX3(i,0,0,numComp,3)] = f_010[i]*tmp0_2 + f_011[i]*tmp0_0 + f_110[i]*tmp0_1 + f_111[i]*tmp0_3;                              const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
480                              o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_11 + f_010[i]*tmp0_9 + f_101[i]*tmp0_13 + f_111[i]*tmp0_10 + tmp0_12*(f_001[i] + f_100[i]) + tmp0_8*(f_011[i] + f_110[i]);                              const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;
481                              o[INDEX3(i,2,0,numComp,3)] = f_010[i]*tmp0_30 + f_011[i]*tmp0_28 + f_110[i]*tmp0_29 + f_111[i]*tmp0_31;                              const double V2=((f_011[i]-f_010[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;
482                              o[INDEX3(i,0,1,numComp,3)] = f_010[i]*tmp0_2 + f_011[i]*tmp0_0 + f_110[i]*tmp0_1 + f_111[i]*tmp0_3;                              const double V3=((f_011[i]-f_010[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;
483                              o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_18 + f_001[i]*tmp0_19 + f_010[i]*tmp0_16 + f_011[i]*tmp0_14 + f_100[i]*tmp0_21 + f_101[i]*tmp0_20 + f_110[i]*tmp0_15 + f_111[i]*tmp0_17;                              o[INDEX3(i,0,0,numComp,3)] = V0;
484                              o[INDEX3(i,2,1,numComp,3)] = f_010[i]*tmp0_34 + f_011[i]*tmp0_32 + f_110[i]*tmp0_33 + f_111[i]*tmp0_35;                              o[INDEX3(i,1,0,numComp,3)] = ((f_010[i]-f_000[i])*C5 + (f_111[i]-f_101[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
485                              o[INDEX3(i,0,2,numComp,3)] = f_010[i]*tmp0_6 + f_011[i]*tmp0_4 + f_110[i]*tmp0_5 + f_111[i]*tmp0_7;                              o[INDEX3(i,2,0,numComp,3)] = V2;
486                              o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_18 + f_001[i]*tmp0_21 + f_010[i]*tmp0_16 + f_011[i]*tmp0_15 + f_100[i]*tmp0_19 + f_101[i]*tmp0_20 + f_110[i]*tmp0_14 + f_111[i]*tmp0_17;                              o[INDEX3(i,0,1,numComp,3)] = V0;
487                              o[INDEX3(i,2,2,numComp,3)] = f_010[i]*tmp0_30 + f_011[i]*tmp0_28 + f_110[i]*tmp0_29 + f_111[i]*tmp0_31;                              o[INDEX3(i,1,1,numComp,3)] = ((f_110[i]-f_100[i])*C5 + (f_011[i]-f_001[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
488                              o[INDEX3(i,0,3,numComp,3)] = f_010[i]*tmp0_6 + f_011[i]*tmp0_4 + f_110[i]*tmp0_5 + f_111[i]*tmp0_7;                              o[INDEX3(i,2,1,numComp,3)] = V3;
489                              o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_25 + f_010[i]*tmp0_23 + f_101[i]*tmp0_27 + f_111[i]*tmp0_24 + tmp0_22*(f_011[i] + f_110[i]) + tmp0_26*(f_001[i] + f_100[i]);                              o[INDEX3(i,0,2,numComp,3)] = V1;
490                              o[INDEX3(i,2,3,numComp,3)] = f_010[i]*tmp0_34 + f_011[i]*tmp0_32 + f_110[i]*tmp0_33 + f_111[i]*tmp0_35;                              o[INDEX3(i,1,2,numComp,3)] = ((f_011[i]-f_001[i])*C5 + (f_110[i]-f_100[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
491                                o[INDEX3(i,2,2,numComp,3)] = V2;
492                                o[INDEX3(i,0,3,numComp,3)] = V1;
493                                o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
494                                o[INDEX3(i,2,3,numComp,3)] = V3;
495                          } /* end of component loop i */                          } /* end of component loop i */
496                      } /* end of k0 loop */                      } /* end of k0 loop */
497                  } /* end of k2 loop */                  } /* end of k2 loop */
498              } /* end of face 3 */              } /* end of face 3 */
499              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
                 const double tmp0_22 = -0.16666666666666666667/h2;  
                 const double tmp0_16 = -0.62200846792814621559/h2;  
                 const double tmp0_33 = 0.044658198738520451079/h2;  
                 const double tmp0_0 = -0.78867513459481288225/h0;  
                 const double tmp0_21 = 0.044658198738520451079/h2;  
                 const double tmp0_17 = -0.16666666666666666667/h2;  
                 const double tmp0_1 = 0.78867513459481288225/h0;  
                 const double tmp0_20 = 0.16666666666666666667/h2;  
                 const double tmp0_14 = 0.78867513459481288225/h1;  
                 const double tmp0_2 = 0.21132486540518711775/h0;  
                 const double tmp0_27 = 0.62200846792814621559/h2;  
                 const double tmp0_15 = 0.21132486540518711775/h1;  
                 const double tmp0_3 = -0.21132486540518711775/h0;  
                 const double tmp0_26 = 0.16666666666666666667/h2;  
                 const double tmp0_12 = -0.21132486540518711775/h1;  
                 const double tmp0_25 = -0.044658198738520451079/h2;  
                 const double tmp0_13 = -0.78867513459481288225/h1;  
                 const double tmp0_24 = -0.16666666666666666667/h2;  
                 const double tmp0_10 = 0.21132486540518711775/h1;  
                 const double tmp0_11 = 0.78867513459481288225/h1;  
                 const double tmp0_34 = 0.16666666666666666667/h2;  
                 const double tmp0_35 = 0.62200846792814621559/h2;  
                 const double tmp0_8 = -0.78867513459481288225/h1;  
                 const double tmp0_29 = 0.16666666666666666667/h2;  
                 const double tmp0_9 = -0.21132486540518711775/h1;  
                 const double tmp0_30 = -0.044658198738520451079/h2;  
                 const double tmp0_28 = 0.044658198738520451079/h2;  
                 const double tmp0_32 = -0.62200846792814621559/h2;  
                 const double tmp0_31 = -0.16666666666666666667/h2;  
                 const double tmp0_18 = -0.044658198738520451079/h2;  
                 const double tmp0_4 = -0.21132486540518711775/h0;  
                 const double tmp0_19 = 0.62200846792814621559/h2;  
                 const double tmp0_5 = 0.21132486540518711775/h0;  
                 const double tmp0_6 = 0.78867513459481288225/h0;  
                 const double tmp0_23 = -0.62200846792814621559/h2;  
                 const double tmp0_7 = -0.78867513459481288225/h0;  
500  #pragma omp for nowait  #pragma omp for nowait
501                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
502                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 721  void Brick::setToGradient(escript::Data& Line 510  void Brick::setToGradient(escript::Data&
510                          const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));                          const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
511                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
512                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
513                              o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_010[i]*tmp0_3 + f_100[i]*tmp0_1 + f_110[i]*tmp0_2;                              const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / h0;
514                              o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_8 + f_010[i]*tmp0_11 + f_100[i]*tmp0_9 + f_110[i]*tmp0_10;                              const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / h0;
515                              o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_16 + f_001[i]*tmp0_19 + f_110[i]*tmp0_18 + f_111[i]*tmp0_21 + tmp0_17*(f_010[i] + f_100[i]) + tmp0_20*(f_011[i] + f_101[i]);                              const double V2=((f_010[i]-f_000[i])*C6 + (f_110[i]-f_100[i])*C2) / h1;
516                              o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_0 + f_010[i]*tmp0_3 + f_100[i]*tmp0_1 + f_110[i]*tmp0_2;                              const double V3=((f_010[i]-f_000[i])*C2 + (f_110[i]-f_100[i])*C6) / h1;
517                              o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_12 + f_010[i]*tmp0_15 + f_100[i]*tmp0_13 + f_110[i]*tmp0_14;                              o[INDEX3(i,0,0,numComp,3)] = V0;
518                              o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_22 + f_001[i]*tmp0_26 + f_010[i]*tmp0_25 + f_011[i]*tmp0_28 + f_100[i]*tmp0_23 + f_101[i]*tmp0_27 + f_110[i]*tmp0_24 + f_111[i]*tmp0_29;                              o[INDEX3(i,1,0,numComp,3)] = V2;
519                              o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_4 + f_010[i]*tmp0_7 + f_100[i]*tmp0_5 + f_110[i]*tmp0_6;                              o[INDEX3(i,2,0,numComp,3)] = ((f_001[i]-f_000[i])*C5 + (f_111[i]-f_110[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
520                              o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_8 + f_010[i]*tmp0_11 + f_100[i]*tmp0_9 + f_110[i]*tmp0_10;                              o[INDEX3(i,0,1,numComp,3)] = V0;
521                              o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_22 + f_001[i]*tmp0_26 + f_010[i]*tmp0_23 + f_011[i]*tmp0_27 + f_100[i]*tmp0_25 + f_101[i]*tmp0_28 + f_110[i]*tmp0_24 + f_111[i]*tmp0_29;                              o[INDEX3(i,1,1,numComp,3)] = V3;
522                              o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_4 + f_010[i]*tmp0_7 + f_100[i]*tmp0_5 + f_110[i]*tmp0_6;                              o[INDEX3(i,2,1,numComp,3)] = ((f_101[i]-f_100[i])*C5 + (f_011[i]-f_010[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
523                              o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_12 + f_010[i]*tmp0_15 + f_100[i]*tmp0_13 + f_110[i]*tmp0_14;                              o[INDEX3(i,0,2,numComp,3)] = V1;
524                              o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_30 + f_001[i]*tmp0_33 + f_110[i]*tmp0_32 + f_111[i]*tmp0_35 + tmp0_31*(f_010[i] + f_100[i]) + tmp0_34*(f_011[i] + f_101[i]);                              o[INDEX3(i,1,2,numComp,3)] = V2;
525                                o[INDEX3(i,2,2,numComp,3)] = ((f_011[i]-f_010[i])*C5 + (f_101[i]-f_100[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
526                                o[INDEX3(i,0,3,numComp,3)] = V1;
527                                o[INDEX3(i,1,3,numComp,3)] = V3;
528                                o[INDEX3(i,2,3,numComp,3)] = ((f_111[i]-f_110[i])*C5 + (f_001[i]-f_000[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
529                          } /* end of component loop i */                          } /* end of component loop i */
530                      } /* end of k0 loop */                      } /* end of k0 loop */
531                  } /* end of k1 loop */                  } /* end of k1 loop */
532              } /* end of face 4 */              } /* end of face 4 */
533              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
                 const double tmp0_22 = 0.16666666666666666667/h2;  
                 const double tmp0_16 = 0.62200846792814621559/h2;  
                 const double tmp0_33 = -0.044658198738520451079/h2;  
                 const double tmp0_0 = -0.78867513459481288225/h0;  
                 const double tmp0_21 = -0.16666666666666666667/h2;  
                 const double tmp0_17 = 0.16666666666666666667/h2;  
                 const double tmp0_1 = 0.78867513459481288225/h0;  
                 const double tmp0_20 = -0.044658198738520451079/h2;  
                 const double tmp0_14 = 0.21132486540518711775/h1;  
                 const double tmp0_2 = -0.21132486540518711775/h0;  
                 const double tmp0_27 = -0.16666666666666666667/h2;  
                 const double tmp0_15 = 0.78867513459481288225/h1;  
                 const double tmp0_3 = 0.21132486540518711775/h0;  
                 const double tmp0_26 = -0.16666666666666666667/h2;  
                 const double tmp0_12 = -0.21132486540518711775/h1;  
                 const double tmp0_25 = 0.16666666666666666667/h2;  
                 const double tmp0_13 = -0.78867513459481288225/h1;  
                 const double tmp0_24 = 0.044658198738520451079/h2;  
                 const double tmp0_10 = 0.78867513459481288225/h1;  
                 const double tmp0_11 = 0.21132486540518711775/h1;  
                 const double tmp0_34 = -0.62200846792814621559/h2;  
                 const double tmp0_35 = -0.16666666666666666667/h2;  
                 const double tmp0_8 = -0.78867513459481288225/h1;  
                 const double tmp0_29 = -0.62200846792814621559/h2;  
                 const double tmp0_9 = -0.21132486540518711775/h1;  
                 const double tmp0_30 = 0.044658198738520451079/h2;  
                 const double tmp0_28 = -0.044658198738520451079/h2;  
                 const double tmp0_32 = 0.62200846792814621559/h2;  
                 const double tmp0_31 = 0.16666666666666666667/h2;  
                 const double tmp0_18 = 0.044658198738520451079/h2;  
                 const double tmp0_4 = -0.21132486540518711775/h0;  
                 const double tmp0_19 = -0.62200846792814621559/h2;  
                 const double tmp0_5 = 0.21132486540518711775/h0;  
                 const double tmp0_6 = -0.78867513459481288225/h0;  
                 const double tmp0_23 = 0.62200846792814621559/h2;  
                 const double tmp0_7 = 0.78867513459481288225/h0;  
534  #pragma omp for nowait  #pragma omp for nowait
535                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
536                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 787  void Brick::setToGradient(escript::Data& Line 544  void Brick::setToGradient(escript::Data&
544                          const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));                          const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
545                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
546                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
547                              o[INDEX3(i,0,0,numComp,3)] = f_001[i]*tmp0_0 + f_011[i]*tmp0_2 + f_101[i]*tmp0_1 + f_111[i]*tmp0_3;                              const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
548                              o[INDEX3(i,1,0,numComp,3)] = f_001[i]*tmp0_8 + f_011[i]*tmp0_10 + f_101[i]*tmp0_9 + f_111[i]*tmp0_11;                              const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;
549                              o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_19 + f_001[i]*tmp0_16 + f_110[i]*tmp0_20 + f_111[i]*tmp0_18 + tmp0_17*(f_011[i] + f_101[i]) + tmp0_21*(f_010[i] + f_100[i]);                              const double V2=((f_011[i]-f_001[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;
550                              o[INDEX3(i,0,1,numComp,3)] = f_001[i]*tmp0_0 + f_011[i]*tmp0_2 + f_101[i]*tmp0_1 + f_111[i]*tmp0_3;                              const double V3=((f_011[i]-f_001[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;
551                              o[INDEX3(i,1,1,numComp,3)] = f_001[i]*tmp0_12 + f_011[i]*tmp0_14 + f_101[i]*tmp0_13 + f_111[i]*tmp0_15;                              o[INDEX3(i,0,0,numComp,3)] = V0;
552                              o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_22 + f_010[i]*tmp0_28 + f_011[i]*tmp0_24 + f_100[i]*tmp0_29 + f_101[i]*tmp0_23 + f_110[i]*tmp0_27 + f_111[i]*tmp0_25;                              o[INDEX3(i,1,0,numComp,3)] = V2;
553                              o[INDEX3(i,0,2,numComp,3)] = f_001[i]*tmp0_4 + f_011[i]*tmp0_6 + f_101[i]*tmp0_5 + f_111[i]*tmp0_7;                              o[INDEX3(i,2,0,numComp,3)] = ((f_001[i]-f_000[i])*C5 + (f_111[i]-f_110[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
554                              o[INDEX3(i,1,2,numComp,3)] = f_001[i]*tmp0_8 + f_011[i]*tmp0_10 + f_101[i]*tmp0_9 + f_111[i]*tmp0_11;                              o[INDEX3(i,0,1,numComp,3)] = V0;
555                              o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_22 + f_010[i]*tmp0_29 + f_011[i]*tmp0_23 + f_100[i]*tmp0_28 + f_101[i]*tmp0_24 + f_110[i]*tmp0_27 + f_111[i]*tmp0_25;                              o[INDEX3(i,1,1,numComp,3)] = V3;
556                              o[INDEX3(i,0,3,numComp,3)] = f_001[i]*tmp0_4 + f_011[i]*tmp0_6 + f_101[i]*tmp0_5 + f_111[i]*tmp0_7;                              o[INDEX3(i,2,1,numComp,3)] = ((f_011[i]-f_010[i])*C0 + (f_101[i]-f_100[i])*C5 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
557                              o[INDEX3(i,1,3,numComp,3)] = f_001[i]*tmp0_12 + f_011[i]*tmp0_14 + f_101[i]*tmp0_13 + f_111[i]*tmp0_15;                              o[INDEX3(i,0,2,numComp,3)] = V1;
558                              o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_33 + f_001[i]*tmp0_30 + f_110[i]*tmp0_34 + f_111[i]*tmp0_32 + tmp0_31*(f_011[i] + f_101[i]) + tmp0_35*(f_010[i] + f_100[i]);                              o[INDEX3(i,1,2,numComp,3)] = V2;
559                                o[INDEX3(i,2,2,numComp,3)] = ((f_011[i]-f_010[i])*C5 + (f_101[i]-f_100[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
560                                o[INDEX3(i,0,3,numComp,3)] = V1;
561                                o[INDEX3(i,1,3,numComp,3)] = V3;
562                                o[INDEX3(i,2,3,numComp,3)] = ((f_001[i]-f_000[i])*C0 + (f_111[i]-f_110[i])*C5 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
563                          } /* end of component loop i */                          } /* end of component loop i */
564                      } /* end of k0 loop */                      } /* end of k0 loop */
565                  } /* end of k1 loop */                  } /* end of k1 loop */
# Line 806  void Brick::setToGradient(escript::Data& Line 567  void Brick::setToGradient(escript::Data&
567          } // end of parallel section          } // end of parallel section
568          /* GENERATOR SNIP_GRAD_FACES BOTTOM */          /* GENERATOR SNIP_GRAD_FACES BOTTOM */
569      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
570          /* GENERATOR SNIP_GRAD_REDUCED_FACES TOP */          /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
571  #pragma omp parallel  #pragma omp parallel
572          {          {
573              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
                 const double tmp0_0 = -0.25/h0;  
                 const double tmp0_4 = -0.5/h2;  
                 const double tmp0_1 = 0.25/h0;  
                 const double tmp0_5 = 0.5/h2;  
                 const double tmp0_2 = -0.5/h1;  
                 const double tmp0_3 = 0.5/h1;  
574  #pragma omp for nowait  #pragma omp for nowait
575                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
576                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 829  void Brick::setToGradient(escript::Data& Line 584  void Brick::setToGradient(escript::Data&
584                          const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));                          const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
585                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
586                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
587                              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]);                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;
588                              o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_001[i]) + tmp0_3*(f_010[i] + f_011[i]);                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / h1;
589                              o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_010[i]) + tmp0_5*(f_001[i] + f_011[i]);                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]-f_000[i]-f_010[i])*C4 / h2;
590                          } /* end of component loop i */                          } /* end of component loop i */
591                      } /* end of k1 loop */                      } /* end of k1 loop */
592                  } /* end of k2 loop */                  } /* end of k2 loop */
593              } /* end of face 0 */              } /* end of face 0 */
594              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
                 const double tmp0_0 = -0.25/h0;  
                 const double tmp0_4 = 0.5/h2;  
                 const double tmp0_1 = 0.25/h0;  
                 const double tmp0_5 = -0.5/h2;  
                 const double tmp0_2 = -0.5/h1;  
                 const double tmp0_3 = 0.5/h1;  
595  #pragma omp for nowait  #pragma omp for nowait
596                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
597                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 856  void Brick::setToGradient(escript::Data& Line 605  void Brick::setToGradient(escript::Data&
605                          const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+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));
606                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
607                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
608                              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]);                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;
609                              o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_100[i] + f_101[i]) + tmp0_3*(f_110[i] + f_111[i]);                              o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / h1;
610                              o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_101[i] + f_111[i]) + tmp0_5*(f_100[i] + f_110[i]);                              o[INDEX3(i,2,0,numComp,3)] = (f_101[i]+f_111[i]-f_100[i]-f_110[i])*C4 / h2;
611                          } /* end of component loop i */                          } /* end of component loop i */
612                      } /* end of k1 loop */                      } /* end of k1 loop */
613                  } /* end of k2 loop */                  } /* end of k2 loop */
614              } /* end of face 1 */              } /* end of face 1 */
615              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
                 const double tmp0_0 = -0.5/h0;  
                 const double tmp0_4 = -0.5/h2;  
                 const double tmp0_1 = 0.5/h0;  
                 const double tmp0_5 = 0.5/h2;  
                 const double tmp0_2 = -0.25/h1;  
                 const double tmp0_3 = 0.25/h1;  
616  #pragma omp for nowait  #pragma omp for nowait
617                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
618                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 883  void Brick::setToGradient(escript::Data& Line 626  void Brick::setToGradient(escript::Data&
626                          const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));                          const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
627                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
628                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
629                              o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_001[i]) + tmp0_1*(f_100[i] + f_101[i]);                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / h0;
630                              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]);                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;
631                              o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_100[i]) + tmp0_5*(f_001[i] + f_101[i]);                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_101[i]-f_000[i]-f_100[i])*C4 / h2;
632                          } /* end of component loop i */                          } /* end of component loop i */
633                      } /* end of k0 loop */                      } /* end of k0 loop */
634                  } /* end of k2 loop */                  } /* end of k2 loop */
635              } /* end of face 2 */              } /* end of face 2 */
636              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
                 const double tmp0_0 = -0.5/h0;  
                 const double tmp0_4 = 0.5/h2;  
                 const double tmp0_1 = 0.5/h0;  
                 const double tmp0_5 = -0.5/h2;  
                 const double tmp0_2 = 0.25/h1;  
                 const double tmp0_3 = -0.25/h1;  
637  #pragma omp for nowait  #pragma omp for nowait
638                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
639                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 910  void Brick::setToGradient(escript::Data& Line 647  void Brick::setToGradient(escript::Data&
647                          const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));                          const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
648                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
649                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
650                              o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_010[i] + f_011[i]) + tmp0_1*(f_110[i] + f_111[i]);                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / h0;
651                              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]);                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;
652                              o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_011[i] + f_111[i]) + tmp0_5*(f_010[i] + f_110[i]);                              o[INDEX3(i,2,0,numComp,3)] = (f_011[i]+f_111[i]-f_010[i]-f_110[i])*C4 / h2;
653                          } /* end of component loop i */                          } /* end of component loop i */
654                      } /* end of k0 loop */                      } /* end of k0 loop */
655                  } /* end of k2 loop */                  } /* end of k2 loop */
656              } /* end of face 3 */              } /* end of face 3 */
657              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
                 const double tmp0_0 = -0.5/h0;  
                 const double tmp0_4 = -0.25/h2;  
                 const double tmp0_1 = 0.5/h0;  
                 const double tmp0_5 = 0.25/h2;  
                 const double tmp0_2 = -0.5/h1;  
                 const double tmp0_3 = 0.5/h1;  
658  #pragma omp for nowait  #pragma omp for nowait
659                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
660                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 937  void Brick::setToGradient(escript::Data& Line 668  void Brick::setToGradient(escript::Data&
668                          const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));                          const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
669                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
670                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
671                              o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_010[i]) + tmp0_1*(f_100[i] + f_110[i]);                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / h0;
672                              o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_100[i]) + tmp0_3*(f_010[i] + f_110[i]);                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / h1;
673                              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]);                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C4 / h2;
674                          } /* end of component loop i */                          } /* end of component loop i */
675                      } /* end of k0 loop */                      } /* end of k0 loop */
676                  } /* end of k1 loop */                  } /* end of k1 loop */
677              } /* end of face 4 */              } /* end of face 4 */
678              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
                 const double tmp0_0 = -0.5/h0;  
                 const double tmp0_4 = 0.25/h2;  
                 const double tmp0_1 = 0.5/h0;  
                 const double tmp0_5 = -0.25/h2;  
                 const double tmp0_2 = -0.5/h1;  
                 const double tmp0_3 = 0.5/h1;  
679  #pragma omp for nowait  #pragma omp for nowait
680                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
681                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 964  void Brick::setToGradient(escript::Data& Line 689  void Brick::setToGradient(escript::Data&
689                          const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));                          const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
690                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
691                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
692                              o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_001[i] + f_011[i]) + tmp0_1*(f_101[i] + f_111[i]);                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / h0;
693                              o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_001[i] + f_101[i]) + tmp0_3*(f_011[i] + f_111[i]);                              o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / h1;
694                              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]);                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / h2;
695                          } /* end of component loop i */                          } /* end of component loop i */
696                      } /* end of k0 loop */                      } /* end of k0 loop */
697                  } /* end of k1 loop */                  } /* end of k1 loop */
# Line 1734  void Brick::interpolateNodesOnElements(e Line 1459  void Brick::interpolateNodesOnElements(e
1459  {  {
1460      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1461      if (reduced) {      if (reduced) {
1462          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */          /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */
1463          const double tmp0_0 = 0.12500000000000000000;          const double c0 = .125;
1464  #pragma omp parallel for  #pragma omp parallel for
1465          for (index_t k2=0; k2 < m_NE2; ++k2) {          for (index_t k2=0; k2 < m_NE2; ++k2) {
1466              for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 1750  void Brick::interpolateNodesOnElements(e Line 1475  void Brick::interpolateNodesOnElements(e
1475                      const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));                      const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1476                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1477                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1478                          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]);                          o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i] + f_100[i] + f_101[i] + f_110[i] + f_111[i]);
1479                      } /* end of component loop i */                      } /* end of component loop i */
1480                  } /* end of k0 loop */                  } /* end of k0 loop */
1481              } /* end of k1 loop */              } /* end of k1 loop */
1482          } /* end of k2 loop */          } /* end of k2 loop */
1483          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */
1484      } else {      } else {
1485          /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */          /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
1486          const double tmp0_3 = 0.0094373878376559314545;          const double c0 = .0094373878376559314545;
1487          const double tmp0_2 = 0.035220810900864519624;          const double c1 = .035220810900864519624;
1488          const double tmp0_1 = 0.13144585576580214704;          const double c2 = .13144585576580214704;
1489          const double tmp0_0 = 0.49056261216234406855;          const double c3 = .49056261216234406855;
1490  #pragma omp parallel for  #pragma omp parallel for
1491          for (index_t k2=0; k2 < m_NE2; ++k2) {          for (index_t k2=0; k2 < m_NE2; ++k2) {
1492              for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 1776  void Brick::interpolateNodesOnElements(e Line 1501  void Brick::interpolateNodesOnElements(e
1501                      const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));                      const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1502                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1503                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1504                          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]);                          o[INDEX2(i,numComp,0)] = f_000[i]*c3 + f_111[i]*c0 + c2*(f_001[i] + f_010[i] + f_100[i]) + c1*(f_011[i] + f_101[i] + f_110[i]);
1505                          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]);                          o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_100[i]*c3 + c2*(f_000[i] + f_101[i] + f_110[i]) + c1*(f_001[i] + f_010[i] + f_111[i]);
1506                          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]);                          o[INDEX2(i,numComp,2)] = f_010[i]*c3 + f_101[i]*c0 + c2*(f_000[i] + f_011[i] + f_110[i]) + c1*(f_001[i] + f_100[i] + f_111[i]);
1507                          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]);                          o[INDEX2(i,numComp,3)] = f_001[i]*c0 + f_110[i]*c3 + c2*(f_010[i] + f_100[i] + f_111[i]) + c1*(f_000[i] + f_011[i] + f_101[i]);
1508                          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]);                          o[INDEX2(i,numComp,4)] = f_001[i]*c3 + f_110[i]*c0 + c2*(f_000[i] + f_011[i] + f_101[i]) + c1*(f_010[i] + f_100[i] + f_111[i]);
1509                          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]);                          o[INDEX2(i,numComp,5)] = f_010[i]*c0 + f_101[i]*c3 + c2*(f_001[i] + f_100[i] + f_111[i]) + c1*(f_000[i] + f_011[i] + f_110[i]);
1510                          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]);                          o[INDEX2(i,numComp,6)] = f_011[i]*c3 + f_100[i]*c0 + c2*(f_001[i] + f_010[i] + f_111[i]) + c1*(f_000[i] + f_101[i] + f_110[i]);
1511                          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]);                          o[INDEX2(i,numComp,7)] = f_000[i]*c0 + f_111[i]*c3 + c2*(f_011[i] + f_101[i] + f_110[i]) + c1*(f_001[i] + f_010[i] + f_100[i]);
1512                      } /* end of component loop i */                      } /* end of component loop i */
1513                  } /* end of k0 loop */                  } /* end of k0 loop */
1514              } /* end of k1 loop */              } /* end of k1 loop */
# Line 1798  void Brick::interpolateNodesOnFaces(escr Line 1523  void Brick::interpolateNodesOnFaces(escr
1523  {  {
1524      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1525      if (reduced) {      if (reduced) {
1526          /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */          const double c0 = .25;
1527          if (m_faceOffset[0] > -1) {  #pragma omp parallel
1528              const double tmp0_0 = 0.25000000000000000000;          {
1529  #pragma omp parallel for              /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */
1530              for (index_t k2=0; k2 < m_NE2; ++k2) {              if (m_faceOffset[0] > -1) {
1531    #pragma omp for nowait
1532                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1533                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1534                            const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
1535                            const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
1536                            const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
1537                            const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
1538                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1539                            for (index_t i=0; i < numComp; ++i) {
1540                                o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]);
1541                            } /* end of component loop i */
1542                        } /* end of k1 loop */
1543                    } /* end of k2 loop */
1544                } /* end of face 0 */
1545                if (m_faceOffset[1] > -1) {
1546    #pragma omp for nowait
1547                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1548                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1549                            const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
1550                            const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
1551                            const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
1552                            const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
1553                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1554                            for (index_t i=0; i < numComp; ++i) {
1555                                o[INDEX2(i,numComp,0)] = c0*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
1556                            } /* end of component loop i */
1557                        } /* end of k1 loop */
1558                    } /* end of k2 loop */
1559                } /* end of face 1 */
1560                if (m_faceOffset[2] > -1) {
1561    #pragma omp for nowait
1562                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1563                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1564                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
1565                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
1566                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
1567                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
1568                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1569                            for (index_t i=0; i < numComp; ++i) {
1570                                o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);
1571                            } /* end of component loop i */
1572                        } /* end of k0 loop */
1573                    } /* end of k2 loop */
1574                } /* end of face 2 */
1575                if (m_faceOffset[3] > -1) {
1576    #pragma omp for nowait
1577                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1578                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1579                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
1580                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
1581                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
1582                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
1583                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1584                            for (index_t i=0; i < numComp; ++i) {
1585                                o[INDEX2(i,numComp,0)] = c0*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);
1586                            } /* end of component loop i */
1587                        } /* end of k0 loop */
1588                    } /* end of k2 loop */
1589                } /* end of face 3 */
1590                if (m_faceOffset[4] > -1) {
1591    #pragma omp for nowait
1592                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1593                      const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1594                      const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));                          const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
1595                      const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));                          const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
1596                      const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));                          const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
1597                      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));
1598                      for (index_t i=0; i < numComp; ++i) {                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1599                          o[INDEX2(i,numComp,0)] = tmp0_0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]);                          for (index_t i=0; i < numComp; ++i) {
1600                      } /* end of component loop i */                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);
1601                            } /* end of component loop i */
1602                        } /* end of k0 loop */
1603                  } /* end of k1 loop */                  } /* end of k1 loop */
1604              } /* end of k2 loop */              } /* end of face 4 */
1605          } /* end of face 0 */              if (m_faceOffset[5] > -1) {
1606          if (m_faceOffset[1] > -1) {  #pragma omp for nowait
             const double tmp0_0 = 0.25000000000000000000;  
 #pragma omp parallel for  
             for (index_t k2=0; k2 < m_NE2; ++k2) {  
1607                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1608                      const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1609                      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));
1610                      const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));                          const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
1611                      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));
1612                      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));
1613                      for (index_t i=0; i < numComp; ++i) {                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1614                          o[INDEX2(i,numComp,0)] = tmp0_0*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);                          for (index_t i=0; i < numComp; ++i) {
1615                      } /* end of component loop i */                              o[INDEX2(i,numComp,0)] = c0*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);
1616                            } /* end of component loop i */
1617                        } /* end of k0 loop */
1618                  } /* end of k1 loop */                  } /* end of k1 loop */
1619              } /* end of k2 loop */              } /* end of face 5 */
1620          } /* end of face 1 */              /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */
1621          if (m_faceOffset[2] > -1) {          } // end of parallel section
             const double tmp0_0 = 0.25000000000000000000;  
 #pragma omp parallel for  
             for (index_t k2=0; k2 < m_NE2; ++k2) {  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));  
                     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));  
                     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));  
                     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX2(i,numComp,0)] = tmp0_0*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k2 loop */  
         } /* end of face 2 */  
         if (m_faceOffset[3] > -1) {  
             const double tmp0_0 = 0.25000000000000000000;  
 #pragma omp parallel for  
             for (index_t k2=0; k2 < m_NE2; ++k2) {  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));  
                     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));  
                     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));  
                     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX2(i,numComp,0)] = tmp0_0*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k2 loop */  
         } /* end of face 3 */  
         if (m_faceOffset[4] > -1) {  
             const double tmp0_0 = 0.25000000000000000000;  
 #pragma omp parallel for  
             for (index_t k1=0; k1 < m_NE1; ++k1) {  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));  
                     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,k1,0, m_N0,m_N1));  
                     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));  
                     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX2(i,numComp,0)] = tmp0_0*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k1 loop */  
         } /* end of face 4 */  
         if (m_faceOffset[5] > -1) {  
             const double tmp0_0 = 0.25000000000000000000;  
 #pragma omp parallel for  
             for (index_t k1=0; k1 < m_NE1; ++k1) {  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));  
                     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));  
                     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));  
                     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));  
                     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX2(i,numComp,0)] = tmp0_0*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k1 loop */  
         } /* end of face 5 */  
         /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */  
1622      } else {      } else {
1623          /* GENERATOR SNIP_INTERPOLATE_FACES TOP */          const double c0 = 0.044658198738520451079;
1624          if (m_faceOffset[0] > -1) {          const double c1 = 0.16666666666666666667;
1625              const double tmp0_2 = 0.044658198738520451079;          const double c2 = 0.62200846792814621559;
1626              const double tmp0_1 = 0.16666666666666666667;  #pragma omp parallel
1627              const double tmp0_0 = 0.62200846792814621559;          {
1628  #pragma omp parallel for              /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */
1629              for (index_t k2=0; k2 < m_NE2; ++k2) {              if (m_faceOffset[0] > -1) {
1630    #pragma omp for nowait
1631                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1632                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1633                            const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
1634                            const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
1635                            const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
1636                            const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
1637                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1638                            for (index_t i=0; i < numComp; ++i) {
1639                                o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_011[i]*c0 + c1*(f_001[i] + f_010[i]);
1640                                o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_010[i]*c2 + c1*(f_000[i] + f_011[i]);
1641                                o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_010[i]*c0 + c1*(f_000[i] + f_011[i]);
1642                                o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_011[i]*c2 + c1*(f_001[i] + f_010[i]);
1643                            } /* end of component loop i */
1644                        } /* end of k1 loop */
1645                    } /* end of k2 loop */
1646                } /* end of face 0 */
1647                if (m_faceOffset[1] > -1) {
1648    #pragma omp for nowait
1649                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1650                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1651                            const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
1652                            const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
1653                            const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
1654                            const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
1655                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1656                            for (index_t i=0; i < numComp; ++i) {
1657                                o[INDEX2(i,numComp,0)] = f_100[i]*c2 + f_111[i]*c0 + c1*(f_101[i] + f_110[i]);
1658                                o[INDEX2(i,numComp,1)] = f_101[i]*c0 + f_110[i]*c2 + c1*(f_100[i] + f_111[i]);
1659                                o[INDEX2(i,numComp,2)] = f_101[i]*c2 + f_110[i]*c0 + c1*(f_100[i] + f_111[i]);
1660                                o[INDEX2(i,numComp,3)] = f_100[i]*c0 + f_111[i]*c2 + c1*(f_101[i] + f_110[i]);
1661                            } /* end of component loop i */
1662                        } /* end of k1 loop */
1663                    } /* end of k2 loop */
1664                } /* end of face 1 */
1665                if (m_faceOffset[2] > -1) {
1666    #pragma omp for nowait
1667                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1668                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1669                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
1670                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
1671                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
1672                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
1673                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1674                            for (index_t i=0; i < numComp; ++i) {
1675                                o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_100[i]);
1676                                o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_101[i]);
1677                                o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_101[i]);
1678                                o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_100[i]);
1679                            } /* end of component loop i */
1680                        } /* end of k0 loop */
1681                    } /* end of k2 loop */
1682                } /* end of face 2 */
1683                if (m_faceOffset[3] > -1) {
1684    #pragma omp for nowait
1685                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1686                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1687                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
1688                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
1689                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
1690                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
1691                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1692                            for (index_t i=0; i < numComp; ++i) {
1693                                o[INDEX2(i,numComp,0)] = f_010[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_110[i]);
1694                                o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_111[i]);
1695                                o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_111[i]);
1696                                o[INDEX2(i,numComp,3)] = f_010[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_110[i]);
1697                            } /* end of component loop i */
1698                        } /* end of k0 loop */
1699                    } /* end of k2 loop */
1700                } /* end of face 3 */
1701                if (m_faceOffset[4] > -1) {
1702    #pragma omp for nowait
1703                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1704                      const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1705                      const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));                          const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
1706                      const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));                          const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
1707                      const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));                          const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
1708                      double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
1709                      for (index_t i=0; i < numComp; ++i) {                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1710                          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) {
1711                          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)] = f_000[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_100[i]);
1712                          o[INDEX2(i,numComp,2)] = f_001[i]*tmp0_0 + f_010[i]*tmp0_2 + tmp0_1*(f_000[i] + f_011[i]);                              o[INDEX2(i,numComp,1)] = f_010[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_110[i]);
1713                          o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_011[i]*tmp0_0 + tmp0_1*(f_001[i] + f_010[i]);                              o[INDEX2(i,numComp,2)] = f_010[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_110[i]);
1714                      } /* end of component loop i */                              o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_100[i]);
1715                            } /* end of component loop i */
1716                        } /* end of k0 loop */
1717                  } /* end of k1 loop */                  } /* end of k1 loop */
1718              } /* end of k2 loop */              } /* end of face 4 */
1719          } /* end of face 0 */              if (m_faceOffset[5] > -1) {
1720          if (m_faceOffset[1] > -1) {  #pragma omp for nowait
             const double tmp0_2 = 0.044658198738520451079;  
             const double tmp0_1 = 0.62200846792814621559;  
             const double tmp0_0 = 0.16666666666666666667;  
 #pragma omp parallel for  
             for (index_t k2=0; k2 < m_NE2; ++k2) {  
1721                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1722                      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) {
1723                      const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));                          const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
1724                      const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));                          const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
1725                      const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));                          const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
1726                      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));
1727                      for (index_t i=0; i < numComp; ++i) {                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1728                          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) {
1729                          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)] = f_001[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_101[i]);
1730                          o[INDEX2(i,numComp,2)] = f_101[i]*tmp0_1 + f_110[i]*tmp0_2 + tmp0_0*(f_100[i] + f_111[i]);                              o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_111[i]);
1731                          o[INDEX2(i,numComp,3)] = f_100[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_101[i] + f_110[i]);                              o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_111[i]);
1732                      } /* end of component loop i */                              o[INDEX2(i,numComp,3)] = f_001[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_101[i]);
1733                            } /* end of component loop i */
1734                        } /* end of k0 loop */
1735                  } /* end of k1 loop */                  } /* end of k1 loop */
1736              } /* end of k2 loop */              } /* end of face 5 */
1737          } /* end of face 1 */              /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
1738          if (m_faceOffset[2] > -1) {          } // end of parallel section
             const double tmp0_2 = 0.044658198738520451079;  
             const double tmp0_1 = 0.16666666666666666667;  
             const double tmp0_0 = 0.62200846792814621559;  
 #pragma omp parallel for  
             for (index_t k2=0; k2 < m_NE2; ++k2) {  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));  
                     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));  
                     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));  
                     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         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,1)] = f_001[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_101[i]);  
                         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,3)] = f_000[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_100[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k2 loop */  
         } /* end of face 2 */  
         if (m_faceOffset[3] > -1) {  
             const double tmp0_2 = 0.044658198738520451079;  
             const double tmp0_1 = 0.62200846792814621559;  
             const double tmp0_0 = 0.16666666666666666667;  
 #pragma omp parallel for  
             for (index_t k2=0; k2 < m_NE2; ++k2) {  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));  
                     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));  
                     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));  
                     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         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,1)] = f_011[i]*tmp0_2 + f_110[i]*tmp0_1 + tmp0_0*(f_010[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]);  
                         o[INDEX2(i,numComp,3)] = f_010[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_011[i] + f_110[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k2 loop */  
         } /* end of face 3 */  
         if (m_faceOffset[4] > -1) {  
             const double tmp0_2 = 0.044658198738520451079;  
             const double tmp0_1 = 0.16666666666666666667;  
             const double tmp0_0 = 0.62200846792814621559;  
 #pragma omp parallel for  
             for (index_t k1=0; k1 < m_NE1; ++k1) {  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));  
                     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));  
                     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));  
                     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));  
                     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         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,1)] = f_010[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_110[i]);  
                         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,3)] = f_000[i]*tmp0_2 + f_110[i]*tmp0_0 + tmp0_1*(f_010[i] + f_100[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k1 loop */  
         } /* end of face 4 */  
         if (m_faceOffset[5] > -1) {  
             const double tmp0_2 = 0.044658198738520451079;  
             const double tmp0_1 = 0.16666666666666666667;  
             const double tmp0_0 = 0.62200846792814621559;  
 #pragma omp parallel for  
             for (index_t k1=0; k1 < m_NE1; ++k1) {  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));  
                     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,k1+1,m_N2-1, m_N0,m_N1));  
                     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));  
                     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         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,1)] = f_011[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_111[i]);  
                         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,3)] = f_001[i]*tmp0_2 + f_111[i]*tmp0_0 + tmp0_1*(f_011[i] + f_101[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k1 loop */  
         } /* end of face 5 */  
         /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */  
1739      }      }
1740  }  }
1741    

Legend:
Removed from v.3722  
changed lines
  Added in v.3731

  ViewVC Help
Powered by ViewVC 1.1.26