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

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

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

revision 4370 by caltinay, Fri Apr 19 06:15:24 2013 UTC revision 4375 by caltinay, Mon Apr 22 05:35:52 2013 UTC
# Line 1095  void Brick::assembleCoordinates(escript: Line 1095  void Brick::assembleCoordinates(escript:
1095  void Brick::assembleGradient(escript::Data& out, escript::Data& in) const  void Brick::assembleGradient(escript::Data& out, escript::Data& in) const
1096  {  {
1097      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
     const double h0 = m_dx[0];  
     const double h1 = m_dx[1];  
     const double h2 = m_dx[2];  
1098      const double C0 = .044658198738520451079;      const double C0 = .044658198738520451079;
1099      const double C1 = .16666666666666666667;      const double C1 = .16666666666666666667;
1100      const double C2 = .21132486540518711775;      const double C2 = .21132486540518711775;
# Line 1132  void Brick::assembleGradient(escript::Da Line 1129  void Brick::assembleGradient(escript::Da
1129                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1130                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));
1131                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1132                              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;                              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) / m_dx[0];
1133                              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;                              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) / m_dx[0];
1134                              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;                              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) / m_dx[0];
1135                              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;                              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) / m_dx[0];
1136                              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;                              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) / m_dx[1];
1137                              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;                              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) / m_dx[1];
1138                              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;                              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) / m_dx[1];
1139                              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;                              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) / m_dx[1];
1140                              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;                              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) / m_dx[2];
1141                              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;                              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) / m_dx[2];
1142                              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;                              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) / m_dx[2];
1143                              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;                              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) / m_dx[2];
1144                              o[INDEX3(i,0,0,numComp,3)] = V0;                              o[INDEX3(i,0,0,numComp,3)] = V0;
1145                              o[INDEX3(i,1,0,numComp,3)] = V4;                              o[INDEX3(i,1,0,numComp,3)] = V4;
1146                              o[INDEX3(i,2,0,numComp,3)] = V8;                              o[INDEX3(i,2,0,numComp,3)] = V8;
# Line 1199  void Brick::assembleGradient(escript::Da Line 1196  void Brick::assembleGradient(escript::Da
1196                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1197                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));
1198                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1199                              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;                              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 / m_dx[0];
1200                              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;                              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 / m_dx[1];
1201                              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;                              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 / m_dx[2];
1202                          } // end of component loop i                          } // end of component loop i
1203                      } // end of k0 loop                      } // end of k0 loop
1204                  } // end of k1 loop                  } // end of k1 loop
# Line 1233  void Brick::assembleGradient(escript::Da Line 1230  void Brick::assembleGradient(escript::Da
1230                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1231                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1232                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1233                              const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / h1;                              const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / m_dx[1];
1234                              const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / h1;                              const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / m_dx[1];
1235                              const double V2=((f_001[i]-f_000[i])*C6 + (f_010[i]-f_011[i])*C2) / h2;                              const double V2=((f_001[i]-f_000[i])*C6 + (f_010[i]-f_011[i])*C2) / m_dx[2];
1236                              const double V3=((f_001[i]-f_000[i])*C2 + (f_011[i]-f_010[i])*C6) / h2;                              const double V3=((f_001[i]-f_000[i])*C2 + (f_011[i]-f_010[i])*C6) / m_dx[2];
1237                              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;                              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) / m_dx[0];
1238                              o[INDEX3(i,1,0,numComp,3)] = V0;                              o[INDEX3(i,1,0,numComp,3)] = V0;
1239                              o[INDEX3(i,2,0,numComp,3)] = V2;                              o[INDEX3(i,2,0,numComp,3)] = V2;
1240                              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;                              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) / m_dx[0];
1241                              o[INDEX3(i,1,1,numComp,3)] = V0;                              o[INDEX3(i,1,1,numComp,3)] = V0;
1242                              o[INDEX3(i,2,1,numComp,3)] = V3;                              o[INDEX3(i,2,1,numComp,3)] = V3;
1243                              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;                              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) / m_dx[0];
1244                              o[INDEX3(i,1,2,numComp,3)] = V1;                              o[INDEX3(i,1,2,numComp,3)] = V1;
1245                              o[INDEX3(i,2,2,numComp,3)] = V2;                              o[INDEX3(i,2,2,numComp,3)] = V2;
1246                              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;                              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) / m_dx[0];
1247                              o[INDEX3(i,1,3,numComp,3)] = V1;                              o[INDEX3(i,1,3,numComp,3)] = V1;
1248                              o[INDEX3(i,2,3,numComp,3)] = V3;                              o[INDEX3(i,2,3,numComp,3)] = V3;
1249                          } // end of component loop i                          } // end of component loop i
# Line 1267  void Brick::assembleGradient(escript::Da Line 1264  void Brick::assembleGradient(escript::Da
1264                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1265                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1266                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1267                              const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;                              const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / m_dx[1];
1268                              const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;                              const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / m_dx[1];
1269                              const double V2=((f_101[i]-f_100[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;                              const double V2=((f_101[i]-f_100[i])*C6 + (f_111[i]-f_110[i])*C2) / m_dx[2];
1270                              const double V3=((f_101[i]-f_100[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;                              const double V3=((f_101[i]-f_100[i])*C2 + (f_111[i]-f_110[i])*C6) / m_dx[2];
1271                              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;                              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) / m_dx[0];
1272                              o[INDEX3(i,1,0,numComp,3)] = V0;                              o[INDEX3(i,1,0,numComp,3)] = V0;
1273                              o[INDEX3(i,2,0,numComp,3)] = V2;                              o[INDEX3(i,2,0,numComp,3)] = V2;
1274                              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;                              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) / m_dx[0];
1275                              o[INDEX3(i,1,1,numComp,3)] = V0;                              o[INDEX3(i,1,1,numComp,3)] = V0;
1276                              o[INDEX3(i,2,1,numComp,3)] = V3;                              o[INDEX3(i,2,1,numComp,3)] = V3;
1277                              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;                              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) / m_dx[0];
1278                              o[INDEX3(i,1,2,numComp,3)] = V1;                              o[INDEX3(i,1,2,numComp,3)] = V1;
1279                              o[INDEX3(i,2,2,numComp,3)] = V2;                              o[INDEX3(i,2,2,numComp,3)] = V2;
1280                              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;                              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) / m_dx[0];
1281                              o[INDEX3(i,1,3,numComp,3)] = V1;                              o[INDEX3(i,1,3,numComp,3)] = V1;
1282                              o[INDEX3(i,2,3,numComp,3)] = V3;                              o[INDEX3(i,2,3,numComp,3)] = V3;
1283                          } // end of component loop i                          } // end of component loop i
# Line 1301  void Brick::assembleGradient(escript::Da Line 1298  void Brick::assembleGradient(escript::Da
1298                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1299                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1300                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1301                              const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / h0;                              const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / m_dx[0];
1302                              const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / h2;                              const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / m_dx[2];
1303                              const double V2=((f_001[i]-f_000[i])*C2 + (f_101[i]-f_100[i])*C6) / h2;                              const double V2=((f_001[i]-f_000[i])*C2 + (f_101[i]-f_100[i])*C6) / m_dx[2];
1304                              o[INDEX3(i,0,0,numComp,3)] = V0;                              o[INDEX3(i,0,0,numComp,3)] = V0;
1305                              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;                              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) / m_dx[1];
1306                              o[INDEX3(i,2,0,numComp,3)] = V1;                              o[INDEX3(i,2,0,numComp,3)] = V1;
1307                              o[INDEX3(i,0,1,numComp,3)] = V0;                              o[INDEX3(i,0,1,numComp,3)] = V0;
1308                              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;                              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) / m_dx[1];
1309                              o[INDEX3(i,2,1,numComp,3)] = V2;                              o[INDEX3(i,2,1,numComp,3)] = V2;
1310                              o[INDEX3(i,0,2,numComp,3)] = V0;                              o[INDEX3(i,0,2,numComp,3)] = V0;
1311                              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;                              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) / m_dx[1];
1312                              o[INDEX3(i,2,2,numComp,3)] = V1;                              o[INDEX3(i,2,2,numComp,3)] = V1;
1313                              o[INDEX3(i,0,3,numComp,3)] = V0;                              o[INDEX3(i,0,3,numComp,3)] = V0;
1314                              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;                              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) / m_dx[1];
1315                              o[INDEX3(i,2,3,numComp,3)] = V2;                              o[INDEX3(i,2,3,numComp,3)] = V2;
1316                          } // end of component loop i                          } // end of component loop i
1317                      } // end of k0 loop                      } // end of k0 loop
# Line 1334  void Brick::assembleGradient(escript::Da Line 1331  void Brick::assembleGradient(escript::Da
1331                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1332                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1333                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1334                              const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;                              const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];
1335                              const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;                              const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];
1336                              const double V2=((f_011[i]-f_010[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;                              const double V2=((f_011[i]-f_010[i])*C6 + (f_111[i]-f_110[i])*C2) / m_dx[2];
1337                              const double V3=((f_011[i]-f_010[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;                              const double V3=((f_011[i]-f_010[i])*C2 + (f_111[i]-f_110[i])*C6) / m_dx[2];
1338                              o[INDEX3(i,0,0,numComp,3)] = V0;                              o[INDEX3(i,0,0,numComp,3)] = V0;
1339                              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;                              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) / m_dx[1];
1340                              o[INDEX3(i,2,0,numComp,3)] = V2;                              o[INDEX3(i,2,0,numComp,3)] = V2;
1341                              o[INDEX3(i,0,1,numComp,3)] = V0;                              o[INDEX3(i,0,1,numComp,3)] = V0;
1342                              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;                              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) / m_dx[1];
1343                              o[INDEX3(i,2,1,numComp,3)] = V3;                              o[INDEX3(i,2,1,numComp,3)] = V3;
1344                              o[INDEX3(i,0,2,numComp,3)] = V1;                              o[INDEX3(i,0,2,numComp,3)] = V1;
1345                              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;                              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) / m_dx[1];
1346                              o[INDEX3(i,2,2,numComp,3)] = V2;                              o[INDEX3(i,2,2,numComp,3)] = V2;
1347                              o[INDEX3(i,0,3,numComp,3)] = V1;                              o[INDEX3(i,0,3,numComp,3)] = V1;
1348                              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;                              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) / m_dx[1];
1349                              o[INDEX3(i,2,3,numComp,3)] = V3;                              o[INDEX3(i,2,3,numComp,3)] = V3;
1350                          } // end of component loop i                          } // end of component loop i
1351                      } // end of k0 loop                      } // end of k0 loop
# Line 1368  void Brick::assembleGradient(escript::Da Line 1365  void Brick::assembleGradient(escript::Da
1365                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1366                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1367                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1368                              const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / h0;                              const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / m_dx[0];
1369                              const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / h0;                              const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / m_dx[0];
1370                              const double V2=((f_010[i]-f_000[i])*C6 + (f_110[i]-f_100[i])*C2) / h1;                              const double V2=((f_010[i]-f_000[i])*C6 + (f_110[i]-f_100[i])*C2) / m_dx[1];
1371                              const double V3=((f_010[i]-f_000[i])*C2 + (f_110[i]-f_100[i])*C6) / h1;                              const double V3=((f_010[i]-f_000[i])*C2 + (f_110[i]-f_100[i])*C6) / m_dx[1];
1372                              o[INDEX3(i,0,0,numComp,3)] = V0;                              o[INDEX3(i,0,0,numComp,3)] = V0;
1373                              o[INDEX3(i,1,0,numComp,3)] = V2;                              o[INDEX3(i,1,0,numComp,3)] = V2;
1374                              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;                              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) / m_dx[2];
1375                              o[INDEX3(i,0,1,numComp,3)] = V0;                              o[INDEX3(i,0,1,numComp,3)] = V0;
1376                              o[INDEX3(i,1,1,numComp,3)] = V3;                              o[INDEX3(i,1,1,numComp,3)] = V3;
1377                              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;                              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) / m_dx[2];
1378                              o[INDEX3(i,0,2,numComp,3)] = V1;                              o[INDEX3(i,0,2,numComp,3)] = V1;
1379                              o[INDEX3(i,1,2,numComp,3)] = V2;                              o[INDEX3(i,1,2,numComp,3)] = V2;
1380                              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;                              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) / m_dx[2];
1381                              o[INDEX3(i,0,3,numComp,3)] = V1;                              o[INDEX3(i,0,3,numComp,3)] = V1;
1382                              o[INDEX3(i,1,3,numComp,3)] = V3;                              o[INDEX3(i,1,3,numComp,3)] = V3;
1383                              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;                              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) / m_dx[2];
1384                          } // end of component loop i                          } // end of component loop i
1385                      } // end of k0 loop                      } // end of k0 loop
1386                  } // end of k1 loop                  } // end of k1 loop
# Line 1402  void Brick::assembleGradient(escript::Da Line 1399  void Brick::assembleGradient(escript::Da
1399                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1400                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1401                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1402                              const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;                              const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];
1403                              const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;                              const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];
1404                              const double V2=((f_011[i]-f_001[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;                              const double V2=((f_011[i]-f_001[i])*C6 + (f_111[i]-f_101[i])*C2) / m_dx[1];
1405                              const double V3=((f_011[i]-f_001[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;                              const double V3=((f_011[i]-f_001[i])*C2 + (f_111[i]-f_101[i])*C6) / m_dx[1];
1406                              o[INDEX3(i,0,0,numComp,3)] = V0;                              o[INDEX3(i,0,0,numComp,3)] = V0;
1407                              o[INDEX3(i,1,0,numComp,3)] = V2;                              o[INDEX3(i,1,0,numComp,3)] = V2;
1408                              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;                              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) / m_dx[2];
1409                              o[INDEX3(i,0,1,numComp,3)] = V0;                              o[INDEX3(i,0,1,numComp,3)] = V0;
1410                              o[INDEX3(i,1,1,numComp,3)] = V3;                              o[INDEX3(i,1,1,numComp,3)] = V3;
1411                              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;                              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) / m_dx[2];
1412                              o[INDEX3(i,0,2,numComp,3)] = V1;                              o[INDEX3(i,0,2,numComp,3)] = V1;
1413                              o[INDEX3(i,1,2,numComp,3)] = V2;                              o[INDEX3(i,1,2,numComp,3)] = V2;
1414                              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;                              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) / m_dx[2];
1415                              o[INDEX3(i,0,3,numComp,3)] = V1;                              o[INDEX3(i,0,3,numComp,3)] = V1;
1416                              o[INDEX3(i,1,3,numComp,3)] = V3;                              o[INDEX3(i,1,3,numComp,3)] = V3;
1417                              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;                              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) / m_dx[2];
1418                          } // end of component loop i                          } // end of component loop i
1419                      } // end of k0 loop                      } // end of k0 loop
1420                  } // end of k1 loop                  } // end of k1 loop
# Line 1449  void Brick::assembleGradient(escript::Da Line 1446  void Brick::assembleGradient(escript::Da
1446                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1447                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1448                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1449                              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;                              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 / m_dx[0];
1450                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / m_dx[1];
1451                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]-f_000[i]-f_010[i])*C4 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]-f_000[i]-f_010[i])*C4 / m_dx[2];
1452                          } // end of component loop i                          } // end of component loop i
1453                      } // end of k1 loop                      } // end of k1 loop
1454                  } // end of k2 loop                  } // end of k2 loop
# Line 1470  void Brick::assembleGradient(escript::Da Line 1467  void Brick::assembleGradient(escript::Da
1467                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1468                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1469                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1470                              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;                              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 / m_dx[0];
1471                              o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / m_dx[1];
1472                              o[INDEX3(i,2,0,numComp,3)] = (f_101[i]+f_111[i]-f_100[i]-f_110[i])*C4 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_101[i]+f_111[i]-f_100[i]-f_110[i])*C4 / m_dx[2];
1473                          } // end of component loop i                          } // end of component loop i
1474                      } // end of k1 loop                      } // end of k1 loop
1475                  } // end of k2 loop                  } // end of k2 loop
# Line 1491  void Brick::assembleGradient(escript::Da Line 1488  void Brick::assembleGradient(escript::Da
1488                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1489                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1490                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1491                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / m_dx[0];
1492                              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;                              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 / m_dx[1];
1493                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_101[i]-f_000[i]-f_100[i])*C4 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_101[i]-f_000[i]-f_100[i])*C4 / m_dx[2];
1494                          } // end of component loop i                          } // end of component loop i
1495                      } // end of k0 loop                      } // end of k0 loop
1496                  } // end of k2 loop                  } // end of k2 loop
# Line 1512  void Brick::assembleGradient(escript::Da Line 1509  void Brick::assembleGradient(escript::Da
1509                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1510                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1511                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1512                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / m_dx[0];
1513                              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;                              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 / m_dx[1];
1514                              o[INDEX3(i,2,0,numComp,3)] = (f_011[i]+f_111[i]-f_010[i]-f_110[i])*C4 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_011[i]+f_111[i]-f_010[i]-f_110[i])*C4 / m_dx[2];
1515                          } // end of component loop i                          } // end of component loop i
1516                      } // end of k0 loop                      } // end of k0 loop
1517                  } // end of k2 loop                  } // end of k2 loop
# Line 1533  void Brick::assembleGradient(escript::Da Line 1530  void Brick::assembleGradient(escript::Da
1530                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1531                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1532                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1533                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / m_dx[0];
1534                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / m_dx[1];
1535                              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;                              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 / m_dx[2];
1536                          } // end of component loop i                          } // end of component loop i
1537                      } // end of k0 loop                      } // end of k0 loop
1538                  } // end of k1 loop                  } // end of k1 loop
# Line 1554  void Brick::assembleGradient(escript::Da Line 1551  void Brick::assembleGradient(escript::Da
1551                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1552                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1553                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1554                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / m_dx[0];
1555                              o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / m_dx[1];
1556                              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;                              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 / m_dx[2];
1557                          } // end of component loop i                          } // end of component loop i
1558                      } // end of k0 loop                      } // end of k0 loop
1559                  } // end of k1 loop                  } // end of k1 loop
# Line 1569  void Brick::assembleGradient(escript::Da Line 1566  void Brick::assembleGradient(escript::Da
1566  void Brick::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const  void Brick::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1567  {  {
1568      const dim_t numComp = arg.getDataPointSize();      const dim_t numComp = arg.getDataPointSize();
     const double h0 = m_dx[0];  
     const double h1 = m_dx[1];  
     const double h2 = m_dx[2];  
1569      const index_t left = (m_offset[0]==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1570      const index_t bottom = (m_offset[1]==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1571      const index_t front = (m_offset[2]==0 ? 0 : 1);      const index_t front = (m_offset[2]==0 ? 0 : 1);
1572      const int fs = arg.getFunctionSpace().getTypeCode();      const int fs = arg.getFunctionSpace().getTypeCode();
1573      if (fs == Elements && arg.actsExpanded()) {      if (fs == Elements && arg.actsExpanded()) {
1574          const double w_0 = h0*h1*h2/8.;          const double w_0 = m_dx[0]*m_dx[1]*m_dx[2]/8.;
1575  #pragma omp parallel  #pragma omp parallel
1576          {          {
1577              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
# Line 1607  void Brick::assembleIntegrate(vector<dou Line 1601  void Brick::assembleIntegrate(vector<dou
1601          } // end of parallel section          } // end of parallel section
1602    
1603      } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {      } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1604          const double w_0 = h0*h1*h2;          const double w_0 = m_dx[0]*m_dx[1]*m_dx[2];
1605  #pragma omp parallel  #pragma omp parallel
1606          {          {
1607              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
# Line 1629  void Brick::assembleIntegrate(vector<dou Line 1623  void Brick::assembleIntegrate(vector<dou
1623          } // end of parallel section          } // end of parallel section
1624    
1625      } else if (fs == FaceElements && arg.actsExpanded()) {      } else if (fs == FaceElements && arg.actsExpanded()) {
1626          const double w_0 = h1*h2/4.;          const double w_0 = m_dx[1]*m_dx[2]/4.;
1627          const double w_1 = h0*h2/4.;          const double w_1 = m_dx[0]*m_dx[2]/4.;
1628          const double w_2 = h0*h1/4.;          const double w_2 = m_dx[0]*m_dx[1]/4.;
1629  #pragma omp parallel  #pragma omp parallel
1630          {          {
1631              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
# Line 1737  void Brick::assembleIntegrate(vector<dou Line 1731  void Brick::assembleIntegrate(vector<dou
1731          } // end of parallel section          } // end of parallel section
1732    
1733      } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {      } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1734          const double w_0 = h1*h2;          const double w_0 = m_dx[1]*m_dx[2];
1735          const double w_1 = h0*h2;          const double w_1 = m_dx[0]*m_dx[2];
1736          const double w_2 = h0*h1;          const double w_2 = m_dx[0]*m_dx[1];
1737  #pragma omp parallel  #pragma omp parallel
1738          {          {
1739              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
# Line 2483  void Brick::interpolateNodesOnElements(e Line 2477  void Brick::interpolateNodesOnElements(e
2477      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
2478      if (reduced) {      if (reduced) {
2479          out.requireWrite();          out.requireWrite();
         const double c0 = .125;  
2480  #pragma omp parallel  #pragma omp parallel
2481          {          {
2482              vector<double> f_000(numComp);              vector<double> f_000(numComp);
# Line 2508  void Brick::interpolateNodesOnElements(e Line 2501  void Brick::interpolateNodesOnElements(e
2501                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2502                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));
2503                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2504                              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]);                              o[INDEX2(i,numComp,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])/8;
2505                          } // end of component loop i                          } // end of component loop i
2506                      } // end of k0 loop                      } // end of k0 loop
2507                  } // end of k1 loop                  } // end of k1 loop
# Line 2567  void Brick::interpolateNodesOnFaces(escr Line 2560  void Brick::interpolateNodesOnFaces(escr
2560      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
2561      if (reduced) {      if (reduced) {
2562          out.requireWrite();          out.requireWrite();
         const double c0 = .25;  
2563  #pragma omp parallel  #pragma omp parallel
2564          {          {
2565              vector<double> f_000(numComp);              vector<double> f_000(numComp);
# Line 2588  void Brick::interpolateNodesOnFaces(escr Line 2580  void Brick::interpolateNodesOnFaces(escr
2580                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2581                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
2582                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2583                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]);                              o[INDEX2(i,numComp,0)] = (f_000[i] + f_001[i] + f_010[i] + f_011[i])/4;
2584                          } // end of component loop i                          } // end of component loop i
2585                      } // end of k1 loop                      } // end of k1 loop
2586                  } // end of k2 loop                  } // end of k2 loop
# Line 2603  void Brick::interpolateNodesOnFaces(escr Line 2595  void Brick::interpolateNodesOnFaces(escr
2595                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2596                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
2597                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2598                              o[INDEX2(i,numComp,0)] = c0*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);                              o[INDEX2(i,numComp,0)] = (f_100[i] + f_101[i] + f_110[i] + f_111[i])/4;
2599                          } // end of component loop i                          } // end of component loop i
2600                      } // end of k1 loop                      } // end of k1 loop
2601                  } // end of k2 loop                  } // end of k2 loop
# Line 2618  void Brick::interpolateNodesOnFaces(escr Line 2610  void Brick::interpolateNodesOnFaces(escr
2610                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2611                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
2612                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2613                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);                              o[INDEX2(i,numComp,0)] = (f_000[i] + f_001[i] + f_100[i] + f_101[i])/4;
2614                          } // end of component loop i                          } // end of component loop i
2615                      } // end of k0 loop                      } // end of k0 loop
2616                  } // end of k2 loop                  } // end of k2 loop
# Line 2633  void Brick::interpolateNodesOnFaces(escr Line 2625  void Brick::interpolateNodesOnFaces(escr
2625                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2626                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
2627                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2628                              o[INDEX2(i,numComp,0)] = c0*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);                              o[INDEX2(i,numComp,0)] = (f_010[i] + f_011[i] + f_110[i] + f_111[i])/4;
2629                          } // end of component loop i                          } // end of component loop i
2630                      } // end of k0 loop                      } // end of k0 loop
2631                  } // end of k2 loop                  } // end of k2 loop
# Line 2648  void Brick::interpolateNodesOnFaces(escr Line 2640  void Brick::interpolateNodesOnFaces(escr
2640                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2641                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
2642                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2643                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);                              o[INDEX2(i,numComp,0)] = (f_000[i] + f_010[i] + f_100[i] + f_110[i])/4;
2644                          } // end of component loop i                          } // end of component loop i
2645                      } // end of k0 loop                      } // end of k0 loop
2646                  } // end of k1 loop                  } // end of k1 loop
# Line 2663  void Brick::interpolateNodesOnFaces(escr Line 2655  void Brick::interpolateNodesOnFaces(escr
2655                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2656                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
2657                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2658                              o[INDEX2(i,numComp,0)] = c0*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);                              o[INDEX2(i,numComp,0)] = (f_001[i] + f_011[i] + f_101[i] + f_111[i])/4;
2659                          } // end of component loop i                          } // end of component loop i
2660                      } // end of k0 loop                      } // end of k0 loop
2661                  } // end of k1 loop                  } // end of k1 loop
# Line 4893  void Brick::assemblePDESingleReduced(Pas Line 4885  void Brick::assemblePDESingleReduced(Pas
4885          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
4886          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
4887  {  {
4888      const double h0 = m_dx[0];      const double w6 = 0.0625*m_dx[0];
4889      const double h1 = m_dx[1];      const double w5 = 0.0625*m_dx[1];
4890      const double h2 = m_dx[2];      const double w1 = 0.0625*m_dx[2];
4891      const double w0 = 0.0625*h1*h2/h0;      const double w14 = 0.03125*m_dx[0]*m_dx[1];
4892      const double w1 = 0.0625*h2;      const double w13 = 0.03125*m_dx[0]*m_dx[2];
4893      const double w2 = -0.0625*h1;      const double w12 = 0.03125*m_dx[1]*m_dx[2];
4894      const double w3 = 0.0625*h0*h2/h1;      const double w18 = 0.015625*m_dx[0]*m_dx[1]*m_dx[2];
4895      const double w4 = -0.0625*h0;      const double w11 = 0.0625*m_dx[0]*m_dx[1]/m_dx[2];
4896      const double w5 = 0.0625*h1;      const double w3 = 0.0625*m_dx[0]*m_dx[2]/m_dx[1];
4897      const double w6 = 0.0625*h0;      const double w0 = 0.0625*m_dx[1]*m_dx[2]/m_dx[0];
     const double w7 = -0.0625*h0*h1/h2;  
     const double w8 = -0.0625*h1*h2/h0;  
     const double w9 = -0.0625*h2;  
     const double w10 = -0.0625*h0*h2/h1;  
     const double w11 = 0.0625*h0*h1/h2;  
     const double w12 = 0.03125*h1*h2;  
     const double w13 = 0.03125*h0*h2;  
     const double w14 = 0.03125*h0*h1;  
     const double w15 = -0.03125*h1*h2;  
     const double w16 = -0.03125*h0*h2;  
     const double w17 = -0.03125*h0*h1;  
     const double w18 = 0.015625*h0*h1*h2;  
     const double w19 = -0.25*h1*h2;  
     const double w20 = -0.25*h0*h2;  
     const double w21 = -0.25*h0*h1;  
     const double w22 = 0.25*h1*h2;  
     const double w23 = 0.25*h0*h2;  
     const double w24 = 0.25*h0*h1;  
     const double w25 = 0.125*h0*h1*h2;  
4898    
4899      rhs.requireWrite();      rhs.requireWrite();
4900  #pragma omp parallel  #pragma omp parallel
# Line 4942  void Brick::assemblePDESingleReduced(Pas Line 4915  void Brick::assemblePDESingleReduced(Pas
4915                          if (!A.isEmpty()) {                          if (!A.isEmpty()) {
4916                              add_EM_S=true;                              add_EM_S=true;
4917                              const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                              const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
4918                              const double A_00 = A_p[INDEX2(0,0,3)];                              const double Aw00 = A_p[INDEX2(0,0,3)]*w0;
4919                              const double A_10 = A_p[INDEX2(1,0,3)];                              const double Aw10 = A_p[INDEX2(1,0,3)]*w1;
4920                              const double A_20 = A_p[INDEX2(2,0,3)];                              const double Aw20 = A_p[INDEX2(2,0,3)]*w5;
4921                              const double A_01 = A_p[INDEX2(0,1,3)];                              const double Aw01 = A_p[INDEX2(0,1,3)]*w1;
4922                              const double A_11 = A_p[INDEX2(1,1,3)];                              const double Aw11 = A_p[INDEX2(1,1,3)]*w3;
4923                              const double A_21 = A_p[INDEX2(2,1,3)];                              const double Aw21 = A_p[INDEX2(2,1,3)]*w6;
4924                              const double A_02 = A_p[INDEX2(0,2,3)];                              const double Aw02 = A_p[INDEX2(0,2,3)]*w5;
4925                              const double A_12 = A_p[INDEX2(1,2,3)];                              const double Aw12 = A_p[INDEX2(1,2,3)]*w6;
4926                              const double A_22 = A_p[INDEX2(2,2,3)];                              const double Aw22 = A_p[INDEX2(2,2,3)]*w11;
4927                              const double tmp0_0 = A_01 + A_10;                              EM_S[INDEX2(0,0,8)]+= Aw00 + Aw01 + Aw02 + Aw10 + Aw11 + Aw12 + Aw20 + Aw21 + Aw22;
4928                              const double tmp1_0 = A_02 + A_20;                              EM_S[INDEX2(1,0,8)]+=-Aw00 - Aw01 - Aw02 + Aw10 + Aw11 + Aw12 + Aw20 + Aw21 + Aw22;
4929                              const double tmp2_0 = A_12 + A_21;                              EM_S[INDEX2(2,0,8)]+= Aw00 + Aw01 + Aw02 - Aw10 - Aw11 - Aw12 + Aw20 + Aw21 + Aw22;
4930                              const double tmp3_1 = A_22*w7;                              EM_S[INDEX2(3,0,8)]+=-Aw00 - Aw01 - Aw02 - Aw10 - Aw11 - Aw12 + Aw20 + Aw21 + Aw22;
4931                              const double tmp10_1 = A_11*w10;                              EM_S[INDEX2(4,0,8)]+= Aw00 + Aw01 + Aw02 + Aw10 + Aw11 + Aw12 - Aw20 - Aw21 - Aw22;
4932                              const double tmp21_1 = A_02*w5;                              EM_S[INDEX2(5,0,8)]+=-Aw00 - Aw01 - Aw02 + Aw10 + Aw11 + Aw12 - Aw20 - Aw21 - Aw22;
4933                              const double tmp2_1 = A_00*w0;                              EM_S[INDEX2(6,0,8)]+= Aw00 + Aw01 + Aw02 - Aw10 - Aw11 - Aw12 - Aw20 - Aw21 - Aw22;
4934                              const double tmp23_1 = tmp2_0*w6;                              EM_S[INDEX2(7,0,8)]+=-Aw00 - Aw01 - Aw02 - Aw10 - Aw11 - Aw12 - Aw20 - Aw21 - Aw22;
4935                              const double tmp19_1 = A_20*w2;                              EM_S[INDEX2(0,1,8)]+=-Aw00 + Aw01 + Aw02 - Aw10 + Aw11 + Aw12 - Aw20 + Aw21 + Aw22;
4936                              const double tmp4_1 = A_11*w3;                              EM_S[INDEX2(1,1,8)]+= Aw00 - Aw01 - Aw02 - Aw10 + Aw11 + Aw12 - Aw20 + Aw21 + Aw22;
4937                              const double tmp22_1 = tmp1_0*w5;                              EM_S[INDEX2(2,1,8)]+=-Aw00 + Aw01 + Aw02 + Aw10 - Aw11 - Aw12 - Aw20 + Aw21 + Aw22;
4938                              const double tmp13_1 = A_21*w4;                              EM_S[INDEX2(3,1,8)]+= Aw00 - Aw01 - Aw02 + Aw10 - Aw11 - Aw12 - Aw20 + Aw21 + Aw22;
4939                              const double tmp5_1 = A_21*w6;                              EM_S[INDEX2(4,1,8)]+=-Aw00 + Aw01 + Aw02 - Aw10 + Aw11 + Aw12 + Aw20 - Aw21 - Aw22;
4940                              const double tmp8_1 = A_00*w8;                              EM_S[INDEX2(5,1,8)]+= Aw00 - Aw01 - Aw02 - Aw10 + Aw11 + Aw12 + Aw20 - Aw21 - Aw22;
4941                              const double tmp7_1 = A_20*w5;                              EM_S[INDEX2(6,1,8)]+=-Aw00 + Aw01 + Aw02 + Aw10 - Aw11 - Aw12 + Aw20 - Aw21 - Aw22;
4942                              const double tmp18_1 = tmp2_0*w4;                              EM_S[INDEX2(7,1,8)]+= Aw00 - Aw01 - Aw02 + Aw10 - Aw11 - Aw12 + Aw20 - Aw21 - Aw22;
4943                              const double tmp6_1 = A_02*w2;                              EM_S[INDEX2(0,2,8)]+= Aw00 - Aw01 + Aw02 + Aw10 - Aw11 + Aw12 + Aw20 - Aw21 + Aw22;
4944                              const double tmp9_1 = A_22*w11;                              EM_S[INDEX2(1,2,8)]+=-Aw00 + Aw01 - Aw02 + Aw10 - Aw11 + Aw12 + Aw20 - Aw21 + Aw22;
4945                              const double tmp15_1 = tmp1_0*w2;                              EM_S[INDEX2(2,2,8)]+= Aw00 - Aw01 + Aw02 - Aw10 + Aw11 - Aw12 + Aw20 - Aw21 + Aw22;
4946                              const double tmp12_1 = A_01*w1;                              EM_S[INDEX2(3,2,8)]+=-Aw00 + Aw01 - Aw02 - Aw10 + Aw11 - Aw12 + Aw20 - Aw21 + Aw22;
4947                              const double tmp0_1 = tmp0_0*w1;                              EM_S[INDEX2(4,2,8)]+= Aw00 - Aw01 + Aw02 + Aw10 - Aw11 + Aw12 - Aw20 + Aw21 - Aw22;
4948                              const double tmp20_1 = A_01*w9;                              EM_S[INDEX2(5,2,8)]+=-Aw00 + Aw01 - Aw02 + Aw10 - Aw11 + Aw12 - Aw20 + Aw21 - Aw22;
4949                              const double tmp14_1 = A_12*w6;                              EM_S[INDEX2(6,2,8)]+= Aw00 - Aw01 + Aw02 - Aw10 + Aw11 - Aw12 - Aw20 + Aw21 - Aw22;
4950                              const double tmp1_1 = A_12*w4;                              EM_S[INDEX2(7,2,8)]+=-Aw00 + Aw01 - Aw02 - Aw10 + Aw11 - Aw12 - Aw20 + Aw21 - Aw22;
4951                              const double tmp16_1 = A_10*w9;                              EM_S[INDEX2(0,3,8)]+=-Aw00 - Aw01 + Aw02 - Aw10 - Aw11 + Aw12 - Aw20 - Aw21 + Aw22;
4952                              const double tmp11_1 = tmp0_0*w9;                              EM_S[INDEX2(1,3,8)]+= Aw00 + Aw01 - Aw02 - Aw10 - Aw11 + Aw12 - Aw20 - Aw21 + Aw22;
4953                              const double tmp17_1 = A_10*w1;                              EM_S[INDEX2(2,3,8)]+=-Aw00 - Aw01 + Aw02 + Aw10 + Aw11 - Aw12 - Aw20 - Aw21 + Aw22;
4954                              EM_S[INDEX2(0,0,8)]+=tmp0_1 + tmp22_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;                              EM_S[INDEX2(3,3,8)]+= Aw00 + Aw01 - Aw02 + Aw10 + Aw11 - Aw12 - Aw20 - Aw21 + Aw22;
4955                              EM_S[INDEX2(1,0,8)]+=tmp17_1 + tmp20_1 + tmp23_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(4,3,8)]+=-Aw00 - Aw01 + Aw02 - Aw10 - Aw11 + Aw12 + Aw20 + Aw21 - Aw22;
4956                              EM_S[INDEX2(2,0,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp1_1 + tmp22_1 + tmp2_1 + tmp5_1 + tmp9_1;                              EM_S[INDEX2(5,3,8)]+= Aw00 + Aw01 - Aw02 - Aw10 - Aw11 + Aw12 + Aw20 + Aw21 - Aw22;
4957                              EM_S[INDEX2(3,0,8)]+=tmp10_1 + tmp11_1 + tmp1_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(6,3,8)]+=-Aw00 - Aw01 + Aw02 + Aw10 + Aw11 - Aw12 + Aw20 + Aw21 - Aw22;
4958                              EM_S[INDEX2(4,0,8)]+=tmp0_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(7,3,8)]+= Aw00 + Aw01 - Aw02 + Aw10 + Aw11 - Aw12 + Aw20 + Aw21 - Aw22;
4959                              EM_S[INDEX2(5,0,8)]+=tmp13_1 + tmp14_1 + tmp15_1 + tmp17_1 + tmp20_1 + tmp3_1 + tmp4_1 + tmp8_1;                              EM_S[INDEX2(0,4,8)]+= Aw00 + Aw01 - Aw02 + Aw10 + Aw11 - Aw12 + Aw20 + Aw21 - Aw22;
4960                              EM_S[INDEX2(6,0,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp18_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(1,4,8)]+=-Aw00 - Aw01 + Aw02 + Aw10 + Aw11 - Aw12 + Aw20 + Aw21 - Aw22;
4961                              EM_S[INDEX2(7,0,8)]+=tmp10_1 + tmp11_1 + tmp15_1 + tmp18_1 + tmp3_1 + tmp8_1;                              EM_S[INDEX2(2,4,8)]+= Aw00 + Aw01 - Aw02 - Aw10 - Aw11 + Aw12 + Aw20 + Aw21 - Aw22;
4962                              EM_S[INDEX2(0,1,8)]+=tmp12_1 + tmp16_1 + tmp19_1 + tmp21_1 + tmp23_1 + tmp4_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(3,4,8)]+=-Aw00 - Aw01 + Aw02 - Aw10 - Aw11 + Aw12 + Aw20 + Aw21 - Aw22;
4963                              EM_S[INDEX2(1,1,8)]+=tmp11_1 + tmp15_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;                              EM_S[INDEX2(4,4,8)]+= Aw00 + Aw01 - Aw02 + Aw10 + Aw11 - Aw12 - Aw20 - Aw21 + Aw22;
4964                              EM_S[INDEX2(2,1,8)]+=tmp0_1 + tmp10_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp5_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(5,4,8)]+=-Aw00 - Aw01 + Aw02 + Aw10 + Aw11 - Aw12 - Aw20 - Aw21 + Aw22;
4965                              EM_S[INDEX2(3,1,8)]+=tmp10_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp20_1 + tmp2_1 + tmp5_1 + tmp9_1;                              EM_S[INDEX2(6,4,8)]+= Aw00 + Aw01 - Aw02 - Aw10 - Aw11 + Aw12 - Aw20 - Aw21 + Aw22;
4966                              EM_S[INDEX2(4,1,8)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp16_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp8_1;                              EM_S[INDEX2(7,4,8)]+=-Aw00 - Aw01 + Aw02 - Aw10 - Aw11 + Aw12 - Aw20 - Aw21 + Aw22;
4967                              EM_S[INDEX2(5,1,8)]+=tmp11_1 + tmp13_1 + tmp14_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(0,5,8)]+=-Aw00 + Aw01 - Aw02 - Aw10 + Aw11 - Aw12 - Aw20 + Aw21 - Aw22;
4968                              EM_S[INDEX2(6,1,8)]+=tmp0_1 + tmp10_1 + tmp18_1 + tmp22_1 + tmp3_1 + tmp8_1;                              EM_S[INDEX2(1,5,8)]+= Aw00 - Aw01 + Aw02 - Aw10 + Aw11 - Aw12 - Aw20 + Aw21 - Aw22;
4969                              EM_S[INDEX2(7,1,8)]+=tmp10_1 + tmp17_1 + tmp18_1 + tmp20_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(2,5,8)]+=-Aw00 + Aw01 - Aw02 + Aw10 - Aw11 + Aw12 - Aw20 + Aw21 - Aw22;
4970                              EM_S[INDEX2(0,2,8)]+=tmp10_1 + tmp13_1 + tmp14_1 + tmp17_1 + tmp20_1 + tmp22_1 + tmp2_1 + tmp9_1;                              EM_S[INDEX2(3,5,8)]+= Aw00 - Aw01 + Aw02 + Aw10 - Aw11 + Aw12 - Aw20 + Aw21 - Aw22;
4971                              EM_S[INDEX2(1,2,8)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp14_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(4,5,8)]+=-Aw00 + Aw01 - Aw02 - Aw10 + Aw11 - Aw12 + Aw20 - Aw21 + Aw22;
4972                              EM_S[INDEX2(2,2,8)]+=tmp11_1 + tmp18_1 + tmp22_1 + tmp2_1 + tmp4_1 + tmp9_1;                              EM_S[INDEX2(5,5,8)]+= Aw00 - Aw01 + Aw02 - Aw10 + Aw11 - Aw12 + Aw20 - Aw21 + Aw22;
4973                              EM_S[INDEX2(3,2,8)]+=tmp12_1 + tmp16_1 + tmp18_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(6,5,8)]+=-Aw00 + Aw01 - Aw02 + Aw10 - Aw11 + Aw12 + Aw20 - Aw21 + Aw22;
4974                              EM_S[INDEX2(4,2,8)]+=tmp10_1 + tmp17_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp23_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(7,5,8)]+= Aw00 - Aw01 + Aw02 + Aw10 - Aw11 + Aw12 + Aw20 - Aw21 + Aw22;
4975                              EM_S[INDEX2(5,2,8)]+=tmp0_1 + tmp10_1 + tmp15_1 + tmp23_1 + tmp3_1 + tmp8_1;                              EM_S[INDEX2(0,6,8)]+= Aw00 - Aw01 - Aw02 + Aw10 - Aw11 - Aw12 + Aw20 - Aw21 - Aw22;
4976                              EM_S[INDEX2(6,2,8)]+=tmp11_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(1,6,8)]+=-Aw00 + Aw01 + Aw02 + Aw10 - Aw11 - Aw12 + Aw20 - Aw21 - Aw22;
4977                              EM_S[INDEX2(7,2,8)]+=tmp12_1 + tmp15_1 + tmp16_1 + tmp1_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;                              EM_S[INDEX2(2,6,8)]+= Aw00 - Aw01 - Aw02 - Aw10 + Aw11 + Aw12 + Aw20 - Aw21 - Aw22;
4978                              EM_S[INDEX2(0,3,8)]+=tmp10_1 + tmp11_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(3,6,8)]+=-Aw00 + Aw01 + Aw02 - Aw10 + Aw11 + Aw12 + Aw20 - Aw21 - Aw22;
4979                              EM_S[INDEX2(1,3,8)]+=tmp10_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp2_1 + tmp9_1;                              EM_S[INDEX2(4,6,8)]+= Aw00 - Aw01 - Aw02 + Aw10 - Aw11 - Aw12 - Aw20 + Aw21 + Aw22;
4980                              EM_S[INDEX2(2,3,8)]+=tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp4_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(5,6,8)]+=-Aw00 + Aw01 + Aw02 + Aw10 - Aw11 - Aw12 - Aw20 + Aw21 + Aw22;
4981                              EM_S[INDEX2(3,3,8)]+=tmp0_1 + tmp15_1 + tmp18_1 + tmp2_1 + tmp4_1 + tmp9_1;                              EM_S[INDEX2(6,6,8)]+= Aw00 - Aw01 - Aw02 - Aw10 + Aw11 + Aw12 - Aw20 + Aw21 + Aw22;
4982                              EM_S[INDEX2(4,3,8)]+=tmp10_1 + tmp11_1 + tmp22_1 + tmp23_1 + tmp3_1 + tmp8_1;                              EM_S[INDEX2(7,6,8)]+=-Aw00 + Aw01 + Aw02 - Aw10 + Aw11 + Aw12 - Aw20 + Aw21 + Aw22;
4983                              EM_S[INDEX2(5,3,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp23_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(0,7,8)]+=-Aw00 - Aw01 - Aw02 - Aw10 - Aw11 - Aw12 - Aw20 - Aw21 - Aw22;
4984                              EM_S[INDEX2(6,3,8)]+=tmp17_1 + tmp1_1 + tmp20_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;                              EM_S[INDEX2(1,7,8)]+= Aw00 + Aw01 + Aw02 - Aw10 - Aw11 - Aw12 - Aw20 - Aw21 - Aw22;
4985                              EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(2,7,8)]+=-Aw00 - Aw01 - Aw02 + Aw10 + Aw11 + Aw12 - Aw20 - Aw21 - Aw22;
4986                              EM_S[INDEX2(0,4,8)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(3,7,8)]+= Aw00 + Aw01 + Aw02 + Aw10 + Aw11 + Aw12 - Aw20 - Aw21 - Aw22;
4987                              EM_S[INDEX2(1,4,8)]+=tmp17_1 + tmp1_1 + tmp20_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;                              EM_S[INDEX2(4,7,8)]+=-Aw00 - Aw01 - Aw02 - Aw10 - Aw11 - Aw12 + Aw20 + Aw21 + Aw22;
4988                              EM_S[INDEX2(2,4,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp23_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(5,7,8)]+= Aw00 + Aw01 + Aw02 - Aw10 - Aw11 - Aw12 + Aw20 + Aw21 + Aw22;
4989                              EM_S[INDEX2(3,4,8)]+=tmp10_1 + tmp11_1 + tmp22_1 + tmp23_1 + tmp3_1 + tmp8_1;                              EM_S[INDEX2(6,7,8)]+=-Aw00 - Aw01 - Aw02 + Aw10 + Aw11 + Aw12 + Aw20 + Aw21 + Aw22;
4990                              EM_S[INDEX2(4,4,8)]+=tmp0_1 + tmp15_1 + tmp18_1 + tmp2_1 + tmp4_1 + tmp9_1;                              EM_S[INDEX2(7,7,8)]+= Aw00 + Aw01 + Aw02 + Aw10 + Aw11 + Aw12 + Aw20 + Aw21 + Aw22;
                             EM_S[INDEX2(5,4,8)]+=tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp4_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(6,4,8)]+=tmp10_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp2_1 + tmp9_1;  
                             EM_S[INDEX2(7,4,8)]+=tmp10_1 + tmp11_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(0,5,8)]+=tmp12_1 + tmp15_1 + tmp16_1 + tmp1_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;  
                             EM_S[INDEX2(1,5,8)]+=tmp11_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(2,5,8)]+=tmp0_1 + tmp10_1 + tmp15_1 + tmp23_1 + tmp3_1 + tmp8_1;  
                             EM_S[INDEX2(3,5,8)]+=tmp10_1 + tmp17_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp23_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(4,5,8)]+=tmp12_1 + tmp16_1 + tmp18_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(5,5,8)]+=tmp11_1 + tmp18_1 + tmp22_1 + tmp2_1 + tmp4_1 + tmp9_1;  
                             EM_S[INDEX2(6,5,8)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp14_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(7,5,8)]+=tmp10_1 + tmp13_1 + tmp14_1 + tmp17_1 + tmp20_1 + tmp22_1 + tmp2_1 + tmp9_1;  
                             EM_S[INDEX2(0,6,8)]+=tmp10_1 + tmp17_1 + tmp18_1 + tmp20_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(1,6,8)]+=tmp0_1 + tmp10_1 + tmp18_1 + tmp22_1 + tmp3_1 + tmp8_1;  
                             EM_S[INDEX2(2,6,8)]+=tmp11_1 + tmp13_1 + tmp14_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(3,6,8)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp16_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp8_1;  
                             EM_S[INDEX2(4,6,8)]+=tmp10_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp20_1 + tmp2_1 + tmp5_1 + tmp9_1;  
                             EM_S[INDEX2(5,6,8)]+=tmp0_1 + tmp10_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp5_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(6,6,8)]+=tmp11_1 + tmp15_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;  
                             EM_S[INDEX2(7,6,8)]+=tmp12_1 + tmp16_1 + tmp19_1 + tmp21_1 + tmp23_1 + tmp4_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(0,7,8)]+=tmp10_1 + tmp11_1 + tmp15_1 + tmp18_1 + tmp3_1 + tmp8_1;  
                             EM_S[INDEX2(1,7,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp18_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(2,7,8)]+=tmp13_1 + tmp14_1 + tmp15_1 + tmp17_1 + tmp20_1 + tmp3_1 + tmp4_1 + tmp8_1;  
                             EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(4,7,8)]+=tmp10_1 + tmp11_1 + tmp1_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(5,7,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp1_1 + tmp22_1 + tmp2_1 + tmp5_1 + tmp9_1;  
                             EM_S[INDEX2(6,7,8)]+=tmp17_1 + tmp20_1 + tmp23_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(7,7,8)]+=tmp0_1 + tmp22_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;  
4991                          }                          }
4992                          ///////////////                          ///////////////
4993                          // process B //                          // process B //
# Line 5049  void Brick::assemblePDESingleReduced(Pas Line 4995  void Brick::assemblePDESingleReduced(Pas
4995                          if (!B.isEmpty()) {                          if (!B.isEmpty()) {
4996                              add_EM_S=true;                              add_EM_S=true;
4997                              const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                              const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
4998                              const double B_0 = B_p[0];                              const double wB0 = B_p[0]*w12;
4999                              const double B_1 = B_p[1];                              const double wB1 = B_p[1]*w13;
5000                              const double B_2 = B_p[2];                              const double wB2 = B_p[2]*w14;
5001                              const double tmp4_1 = B_0*w15;                              EM_S[INDEX2(0,0,8)]+=-wB0 - wB1 - wB2;
5002                              const double tmp3_1 = B_1*w16;                              EM_S[INDEX2(0,1,8)]+=-wB0 - wB1 - wB2;
5003                              const double tmp2_1 = B_0*w12;                              EM_S[INDEX2(0,2,8)]+=-wB0 - wB1 - wB2;
5004                              const double tmp5_1 = B_2*w17;                              EM_S[INDEX2(0,3,8)]+=-wB0 - wB1 - wB2;
5005                              const double tmp1_1 = B_2*w14;                              EM_S[INDEX2(0,4,8)]+=-wB0 - wB1 - wB2;
5006                              const double tmp0_1 = B_1*w13;                              EM_S[INDEX2(0,5,8)]+=-wB0 - wB1 - wB2;
5007                              EM_S[INDEX2(0,0,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,6,8)]+=-wB0 - wB1 - wB2;
5008                              EM_S[INDEX2(1,0,8)]+=tmp2_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(0,7,8)]+=-wB0 - wB1 - wB2;
5009                              EM_S[INDEX2(2,0,8)]+=tmp0_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(1,0,8)]+= wB0 - wB1 - wB2;
5010                              EM_S[INDEX2(3,0,8)]+=tmp0_1 + tmp2_1 + tmp5_1;                              EM_S[INDEX2(1,1,8)]+= wB0 - wB1 - wB2;
5011                              EM_S[INDEX2(4,0,8)]+=tmp1_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(1,2,8)]+= wB0 - wB1 - wB2;
5012                              EM_S[INDEX2(5,0,8)]+=tmp1_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(1,3,8)]+= wB0 - wB1 - wB2;
5013                              EM_S[INDEX2(6,0,8)]+=tmp0_1 + tmp1_1 + tmp4_1;                              EM_S[INDEX2(1,4,8)]+= wB0 - wB1 - wB2;
5014                              EM_S[INDEX2(7,0,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(1,5,8)]+= wB0 - wB1 - wB2;
5015                              EM_S[INDEX2(0,1,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(1,6,8)]+= wB0 - wB1 - wB2;
5016                              EM_S[INDEX2(1,1,8)]+=tmp2_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(1,7,8)]+= wB0 - wB1 - wB2;
5017                              EM_S[INDEX2(2,1,8)]+=tmp0_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(2,0,8)]+=-wB0 + wB1 - wB2;
5018                              EM_S[INDEX2(3,1,8)]+=tmp0_1 + tmp2_1 + tmp5_1;                              EM_S[INDEX2(2,1,8)]+=-wB0 + wB1 - wB2;
5019                              EM_S[INDEX2(4,1,8)]+=tmp1_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(2,2,8)]+=-wB0 + wB1 - wB2;
5020                              EM_S[INDEX2(5,1,8)]+=tmp1_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(2,3,8)]+=-wB0 + wB1 - wB2;
5021                              EM_S[INDEX2(6,1,8)]+=tmp0_1 + tmp1_1 + tmp4_1;                              EM_S[INDEX2(2,4,8)]+=-wB0 + wB1 - wB2;
5022                              EM_S[INDEX2(7,1,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(2,5,8)]+=-wB0 + wB1 - wB2;
5023                              EM_S[INDEX2(0,2,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(2,6,8)]+=-wB0 + wB1 - wB2;
5024                              EM_S[INDEX2(1,2,8)]+=tmp2_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(2,7,8)]+=-wB0 + wB1 - wB2;
5025                              EM_S[INDEX2(2,2,8)]+=tmp0_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(3,0,8)]+= wB0 + wB1 - wB2;
5026                              EM_S[INDEX2(3,2,8)]+=tmp0_1 + tmp2_1 + tmp5_1;                              EM_S[INDEX2(3,1,8)]+= wB0 + wB1 - wB2;
5027                              EM_S[INDEX2(4,2,8)]+=tmp1_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(3,2,8)]+= wB0 + wB1 - wB2;
5028                              EM_S[INDEX2(5,2,8)]+=tmp1_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(3,3,8)]+= wB0 + wB1 - wB2;
5029                              EM_S[INDEX2(6,2,8)]+=tmp0_1 + tmp1_1 + tmp4_1;                              EM_S[INDEX2(3,4,8)]+= wB0 + wB1 - wB2;
5030                              EM_S[INDEX2(7,2,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(3,5,8)]+= wB0 + wB1 - wB2;
5031                              EM_S[INDEX2(0,3,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(3,6,8)]+= wB0 + wB1 - wB2;
5032                              EM_S[INDEX2(1,3,8)]+=tmp2_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(3,7,8)]+= wB0 + wB1 - wB2;
5033                              EM_S[INDEX2(2,3,8)]+=tmp0_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(4,0,8)]+=-wB0 - wB1 + wB2;
5034                              EM_S[INDEX2(3,3,8)]+=tmp0_1 + tmp2_1 + tmp5_1;                              EM_S[INDEX2(4,1,8)]+=-wB0 - wB1 + wB2;
5035                              EM_S[INDEX2(4,3,8)]+=tmp1_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(4,2,8)]+=-wB0 - wB1 + wB2;
5036                              EM_S[INDEX2(5,3,8)]+=tmp1_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(4,3,8)]+=-wB0 - wB1 + wB2;
5037                              EM_S[INDEX2(6,3,8)]+=tmp0_1 + tmp1_1 + tmp4_1;                              EM_S[INDEX2(4,4,8)]+=-wB0 - wB1 + wB2;
5038                              EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(4,5,8)]+=-wB0 - wB1 + wB2;
5039                              EM_S[INDEX2(0,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(4,6,8)]+=-wB0 - wB1 + wB2;
5040                              EM_S[INDEX2(1,4,8)]+=tmp2_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(4,7,8)]+=-wB0 - wB1 + wB2;
5041                              EM_S[INDEX2(2,4,8)]+=tmp0_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(5,0,8)]+= wB0 - wB1 + wB2;
5042                              EM_S[INDEX2(3,4,8)]+=tmp0_1 + tmp2_1 + tmp5_1;                              EM_S[INDEX2(5,1,8)]+= wB0 - wB1 + wB2;
5043                              EM_S[INDEX2(4,4,8)]+=tmp1_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(5,2,8)]+= wB0 - wB1 + wB2;
5044                              EM_S[INDEX2(5,4,8)]+=tmp1_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(5,3,8)]+= wB0 - wB1 + wB2;
5045                              EM_S[INDEX2(6,4,8)]+=tmp0_1 + tmp1_1 + tmp4_1;                              EM_S[INDEX2(5,4,8)]+= wB0 - wB1 + wB2;
5046                              EM_S[INDEX2(7,4,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(5,5,8)]+= wB0 - wB1 + wB2;
5047                              EM_S[INDEX2(0,5,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(5,6,8)]+= wB0 - wB1 + wB2;
5048                              EM_S[INDEX2(1,5,8)]+=tmp2_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(5,7,8)]+= wB0 - wB1 + wB2;
5049                              EM_S[INDEX2(2,5,8)]+=tmp0_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(6,0,8)]+=-wB0 + wB1 + wB2;
5050                              EM_S[INDEX2(3,5,8)]+=tmp0_1 + tmp2_1 + tmp5_1;                              EM_S[INDEX2(6,1,8)]+=-wB0 + wB1 + wB2;
5051                              EM_S[INDEX2(4,5,8)]+=tmp1_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(6,2,8)]+=-wB0 + wB1 + wB2;
5052                              EM_S[INDEX2(5,5,8)]+=tmp1_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(6,3,8)]+=-wB0 + wB1 + wB2;
5053                              EM_S[INDEX2(6,5,8)]+=tmp0_1 + tmp1_1 + tmp4_1;                              EM_S[INDEX2(6,4,8)]+=-wB0 + wB1 + wB2;
5054                              EM_S[INDEX2(7,5,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(6,5,8)]+=-wB0 + wB1 + wB2;
5055                              EM_S[INDEX2(0,6,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(6,6,8)]+=-wB0 + wB1 + wB2;
5056                              EM_S[INDEX2(1,6,8)]+=tmp2_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(6,7,8)]+=-wB0 + wB1 + wB2;
5057                              EM_S[INDEX2(2,6,8)]+=tmp0_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(7,0,8)]+= wB0 + wB1 + wB2;
5058                              EM_S[INDEX2(3,6,8)]+=tmp0_1 + tmp2_1 + tmp5_1;                              EM_S[INDEX2(7,1,8)]+= wB0 + wB1 + wB2;
5059                              EM_S[INDEX2(4,6,8)]+=tmp1_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(7,2,8)]+= wB0 + wB1 + wB2;
5060                              EM_S[INDEX2(5,6,8)]+=tmp1_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(7,3,8)]+= wB0 + wB1 + wB2;
5061                              EM_S[INDEX2(6,6,8)]+=tmp0_1 + tmp1_1 + tmp4_1;                              EM_S[INDEX2(7,4,8)]+= wB0 + wB1 + wB2;
5062                              EM_S[INDEX2(7,6,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(7,5,8)]+= wB0 + wB1 + wB2;
5063                              EM_S[INDEX2(0,7,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(7,6,8)]+= wB0 + wB1 + wB2;
5064                              EM_S[INDEX2(1,7,8)]+=tmp2_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(7,7,8)]+= wB0 + wB1 + wB2;
                             EM_S[INDEX2(2,7,8)]+=tmp0_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp2_1 + tmp5_1;  
                             EM_S[INDEX2(4,7,8)]+=tmp1_1 + tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(5,7,8)]+=tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(6,7,8)]+=tmp0_1 + tmp1_1 + tmp4_1;  
                             EM_S[INDEX2(7,7,8)]+=tmp0_1 + tmp1_1 + tmp2_1;  
5065                          }                          }
5066                          ///////////////                          ///////////////
5067                          // process C //                          // process C //
# Line 5129  void Brick::assemblePDESingleReduced(Pas Line 5069  void Brick::assemblePDESingleReduced(Pas
5069                          if (!C.isEmpty()) {                          if (!C.isEmpty()) {
5070                              add_EM_S=true;                              add_EM_S=true;
5071                              const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                              const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
5072                              const double C_0 = C_p[0];                              const double wC0 = C_p[0]*w12;
5073                              const double C_1 = C_p[1];                              const double wC1 = C_p[1]*w13;
5074                              const double C_2 = C_p[2];                              const double wC2 = C_p[2]*w14;
5075                              const double tmp5_1 = C_0*w15;                              EM_S[INDEX2(0,0,8)]+=-wC0 - wC1 - wC2;
5076                              const double tmp2_1 = C_0*w12;                              EM_S[INDEX2(1,0,8)]+=-wC0 - wC1 - wC2;
5077                              const double tmp4_1 = C_1*w16;                              EM_S[INDEX2(2,0,8)]+=-wC0 - wC1 - wC2;
5078                              const double tmp1_1 = C_2*w17;                              EM_S[INDEX2(3,0,8)]+=-wC0 - wC1 - wC2;
5079                              const double tmp3_1 = C_2*w14;                              EM_S[INDEX2(4,0,8)]+=-wC0 - wC1 - wC2;
5080                              const double tmp0_1 = C_1*w13;                              EM_S[INDEX2(5,0,8)]+=-wC0 - wC1 - wC2;
5081                              EM_S[INDEX2(0,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(6,0,8)]+=-wC0 - wC1 - wC2;
5082                              EM_S[INDEX2(1,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(7,0,8)]+=-wC0 - wC1 - wC2;
5083                              EM_S[INDEX2(2,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,1,8)]+= wC0 - wC1 - wC2;
5084                              EM_S[INDEX2(3,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(1,1,8)]+= wC0 - wC1 - wC2;
5085                              EM_S[INDEX2(4,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(2,1,8)]+= wC0 - wC1 - wC2;
5086                              EM_S[INDEX2(5,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(3,1,8)]+= wC0 - wC1 - wC2;
5087                              EM_S[INDEX2(6,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(4,1,8)]+= wC0 - wC1 - wC2;
5088                              EM_S[INDEX2(7,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(5,1,8)]+= wC0 - wC1 - wC2;
5089                              EM_S[INDEX2(0,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                              EM_S[INDEX2(6,1,8)]+= wC0 - wC1 - wC2;
5090                              EM_S[INDEX2(1,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                              EM_S[INDEX2(7,1,8)]+= wC0 - wC1 - wC2;
5091                              EM_S[INDEX2(2,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                              EM_S[INDEX2(0,2,8)]+=-wC0 + wC1 - wC2;
5092                              EM_S[INDEX2(3,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                              EM_S[INDEX2(1,2,8)]+=-wC0 + wC1 - wC2;
5093                              EM_S[INDEX2(4,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                              EM_S[INDEX2(2,2,8)]+=-wC0 + wC1 - wC2;
5094                              EM_S[INDEX2(5,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                              EM_S[INDEX2(3,2,8)]+=-wC0 + wC1 - wC2;
5095                              EM_S[INDEX2(6,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                              EM_S[INDEX2(4,2,8)]+=-wC0 + wC1 - wC2;
5096                              EM_S[INDEX2(7,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                              EM_S[INDEX2(5,2,8)]+=-wC0 + wC1 - wC2;
5097                              EM_S[INDEX2(0,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                              EM_S[INDEX2(6,2,8)]+=-wC0 + wC1 - wC2;
5098                              EM_S[INDEX2(1,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                              EM_S[INDEX2(7,2,8)]+=-wC0 + wC1 - wC2;
5099                              EM_S[INDEX2(2,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                              EM_S[INDEX2(0,3,8)]+= wC0 + wC1 - wC2;
5100                              EM_S[INDEX2(3,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                              EM_S[INDEX2(1,3,8)]+= wC0 + wC1 - wC2;
5101                              EM_S[INDEX2(4,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                              EM_S[INDEX2(2,3,8)]+= wC0 + wC1 - wC2;
5102                              EM_S[INDEX2(5,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                              EM_S[INDEX2(3,3,8)]+= wC0 + wC1 - wC2;
5103                              EM_S[INDEX2(6,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                              EM_S[INDEX2(4,3,8)]+= wC0 + wC1 - wC2;
5104                              EM_S[INDEX2(7,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                              EM_S[INDEX2(5,3,8)]+= wC0 + wC1 - wC2;
5105                              EM_S[INDEX2(0,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(6,3,8)]+= wC0 + wC1 - wC2;
5106                              EM_S[INDEX2(1,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(7,3,8)]+= wC0 + wC1 - wC2;
5107                              EM_S[INDEX2(2,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(0,4,8)]+=-wC0 - wC1 + wC2;
5108                              EM_S[INDEX2(3,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(1,4,8)]+=-wC0 - wC1 + wC2;
5109                              EM_S[INDEX2(4,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(2,4,8)]+=-wC0 - wC1 + wC2;
5110                              EM_S[INDEX2(5,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(3,4,8)]+=-wC0 - wC1 + wC2;
5111                              EM_S[INDEX2(6,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(4,4,8)]+=-wC0 - wC1 + wC2;
5112                              EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(5,4,8)]+=-wC0 - wC1 + wC2;
5113                              EM_S[INDEX2(0,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(6,4,8)]+=-wC0 - wC1 + wC2;
5114                              EM_S[INDEX2(1,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(7,4,8)]+=-wC0 - wC1 + wC2;
5115                              EM_S[INDEX2(2,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,5,8)]+= wC0 - wC1 + wC2;
5116                              EM_S[INDEX2(3,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(1,5,8)]+= wC0 - wC1 + wC2;
5117                              EM_S[INDEX2(4,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(2,5,8)]+= wC0 - wC1 + wC2;
5118                              EM_S[INDEX2(5,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(3,5,8)]+= wC0 - wC1 + wC2;
5119                              EM_S[INDEX2(6,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(4,5,8)]+= wC0 - wC1 + wC2;
5120                              EM_S[INDEX2(7,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_S[INDEX2(5,5,8)]+= wC0 - wC1 + wC2;
5121                              EM_S[INDEX2(0,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(6,5,8)]+= wC0 - wC1 + wC2;
5122                              EM_S[INDEX2(1,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(7,5,8)]+= wC0 - wC1 + wC2;
5123                              EM_S[INDEX2(2,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(0,6,8)]+=-wC0 + wC1 + wC2;
5124                              EM_S[INDEX2(3,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(1,6,8)]+=-wC0 + wC1 + wC2;
5125                              EM_S[INDEX2(4,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(2,6,8)]+=-wC0 + wC1 + wC2;
5126                              EM_S[INDEX2(5,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(3,6,8)]+=-wC0 + wC1 + wC2;
5127                              EM_S[INDEX2(6,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(4,6,8)]+=-wC0 + wC1 + wC2;
5128                              EM_S[INDEX2(7,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                              EM_S[INDEX2(5,6,8)]+=-wC0 + wC1 + wC2;
5129                              EM_S[INDEX2(0,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(6,6,8)]+=-wC0 + wC1 + wC2;
5130                              EM_S[INDEX2(1,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(7,6,8)]+=-wC0 + wC1 + wC2;
5131                              EM_S[INDEX2(2,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(0,7,8)]+= wC0 + wC1 + wC2;
5132                              EM_S[INDEX2(3,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(1,7,8)]+= wC0 + wC1 + wC2;
5133                              EM_S[INDEX2(4,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(2,7,8)]+= wC0 + wC1 + wC2;
5134                              EM_S[INDEX2(5,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(3,7,8)]+= wC0 + wC1 + wC2;
5135                              EM_S[INDEX2(6,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(4,7,8)]+= wC0 + wC1 + wC2;
5136                              EM_S[INDEX2(7,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                              EM_S[INDEX2(5,7,8)]+= wC0 + wC1 + wC2;
5137                              EM_S[INDEX2(0,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(6,7,8)]+= wC0 + wC1 + wC2;
5138                              EM_S[INDEX2(1,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(7,7,8)]+= wC0 + wC1 + wC2;
                             EM_S[INDEX2(2,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(4,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(5,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(6,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(7,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;  
5139                          }                          }
5140                          ///////////////                          ///////////////
5141                          // process D //                          // process D //
# Line 5209  void Brick::assemblePDESingleReduced(Pas Line 5143  void Brick::assemblePDESingleReduced(Pas
5143                          if (!D.isEmpty()) {                          if (!D.isEmpty()) {
5144                              add_EM_S=true;                              add_EM_S=true;
5145                              const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                              const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
5146                              const double tmp0_1 = D_p[0]*w18;                              EM_S[INDEX2(0,0,8)]+=D_p[0]*w18;
5147                              EM_S[INDEX2(0,0,8)]+=tmp0_1;                              EM_S[INDEX2(1,0,8)]+=D_p[0]*w18;
5148                              EM_S[INDEX2(1,0,8)]+=tmp0_1;                              EM_S[INDEX2(2,0,8)]+=D_p[0]*w18;
5149                              EM_S[INDEX2(2,0,8)]+=tmp0_1;                              EM_S[INDEX2(3,0,8)]+=D_p[0]*w18;
5150                              EM_S[INDEX2(3,0,8)]+=tmp0_1;                              EM_S[INDEX2(4,0,8)]+=D_p[0]*w18;
5151                              EM_S[INDEX2(4,0,8)]+=tmp0_1;                              EM_S[INDEX2(5,0,8)]+=D_p[0]*w18;
5152                              EM_S[INDEX2(5,0,8)]+=tmp0_1;                              EM_S[INDEX2(6,0,8)]+=D_p[0]*w18;
5153                              EM_S[INDEX2(6,0,8)]+=tmp0_1;                              EM_S[INDEX2(7,0,8)]+=D_p[0]*w18;
5154                              EM_S[INDEX2(7,0,8)]+=tmp0_1;                              EM_S[INDEX2(0,1,8)]+=D_p[0]*w18;
5155                              EM_S[INDEX2(0,1,8)]+=tmp0_1;                              EM_S[INDEX2(1,1,8)]+=D_p[0]*w18;
5156                              EM_S[INDEX2(1,1,8)]+=tmp0_1;                              EM_S[INDEX2(2,1,8)]+=D_p[0]*w18;
5157                              EM_S[INDEX2(2,1,8)]+=tmp0_1;                              EM_S[INDEX2(3,1,8)]+=D_p[0]*w18;
5158                              EM_S[INDEX2(3,1,8)]+=tmp0_1;                              EM_S[INDEX2(4,1,8)]+=D_p[0]*w18;
5159                              EM_S[INDEX2(4,1,8)]+=tmp0_1;                              EM_S[INDEX2(5,1,8)]+=D_p[0]*w18;
5160                              EM_S[INDEX2(5,1,8)]+=tmp0_1;                              EM_S[INDEX2(6,1,8)]+=D_p[0]*w18;
5161                              EM_S[INDEX2(6,1,8)]+=tmp0_1;                              EM_S[INDEX2(7,1,8)]+=D_p[0]*w18;
5162                              EM_S[INDEX2(7,1,8)]+=tmp0_1;                              EM_S[INDEX2(0,2,8)]+=D_p[0]*w18;
5163                              EM_S[INDEX2(0,2,8)]+=tmp0_1;                              EM_S[INDEX2(1,2,8)]+=D_p[0]*w18;
5164                              EM_S[INDEX2(1,2,8)]+=tmp0_1;                              EM_S[INDEX2(2,2,8)]+=D_p[0]*w18;
5165                              EM_S[INDEX2(2,2,8)]+=tmp0_1;                              EM_S[INDEX2(3,2,8)]+=D_p[0]*w18;
5166                              EM_S[INDEX2(3,2,8)]+=tmp0_1;                              EM_S[INDEX2(4,2,8)]+=D_p[0]*w18;
5167                              EM_S[INDEX2(4,2,8)]+=tmp0_1;                              EM_S[INDEX2(5,2,8)]+=D_p[0]*w18;
5168                              EM_S[INDEX2(5,2,8)]+=tmp0_1;                              EM_S[INDEX2(6,2,8)]+=D_p[0]*w18;
5169                              EM_S[INDEX2(6,2,8)]+=tmp0_1;                              EM_S[INDEX2(7,2,8)]+=D_p[0]*w18;
5170                              EM_S[INDEX2(7,2,8)]+=tmp0_1;                              EM_S[INDEX2(0,3,8)]+=D_p[0]*w18;
5171                              EM_S[INDEX2(0,3,8)]+=tmp0_1;                              EM_S[INDEX2(1,3,8)]+=D_p[0]*w18;
5172                              EM_S[INDEX2(1,3,8)]+=tmp0_1;                              EM_S[INDEX2(2,3,8)]+=D_p[0]*w18;
5173                              EM_S[INDEX2(2,3,8)]+=tmp0_1;                              EM_S[INDEX2(3,3,8)]+=D_p[0]*w18;
5174                              EM_S[INDEX2(3,3,8)]+=tmp0_1;                              EM_S[INDEX2(4,3,8)]+=D_p[0]*w18;
5175                              EM_S[INDEX2(4,3,8)]+=tmp0_1;                              EM_S[INDEX2(5,3,8)]+=D_p[0]*w18;
5176                              EM_S[INDEX2(5,3,8)]+=tmp0_1;                              EM_S[INDEX2(6,3,8)]+=D_p[0]*w18;
5177                              EM_S[INDEX2(6,3,8)]+=tmp0_1;                              EM_S[INDEX2(7,3,8)]+=D_p[0]*w18;
5178                              EM_S[INDEX2(7,3,8)]+=tmp0_1;                              EM_S[INDEX2(0,4,8)]+=D_p[0]*w18;
5179                              EM_S[INDEX2(0,4,8)]+=tmp0_1;                              EM_S[INDEX2(1,4,8)]+=D_p[0]*w18;
5180                              EM_S[INDEX2(1,4,8)]+=tmp0_1;                              EM_S[INDEX2(2,4,8)]+=D_p[0]*w18;
5181                              EM_S[INDEX2(2,4,8)]+=tmp0_1;                              EM_S[INDEX2(3,4,8)]+=D_p[0]*w18;
5182                              EM_S[INDEX2(3,4,8)]+=tmp0_1;                              EM_S[INDEX2(4,4,8)]+=D_p[0]*w18;
5183                              EM_S[INDEX2(4,4,8)]+=tmp0_1;                              EM_S[INDEX2(5,4,8)]+=D_p[0]*w18;
5184                              EM_S[INDEX2(5,4,8)]+=tmp0_1;                              EM_S[INDEX2(6,4,8)]+=D_p[0]*w18;
5185                              EM_S[INDEX2(6,4,8)]+=tmp0_1;                              EM_S[INDEX2(7,4,8)]+=D_p[0]*w18;
5186                              EM_S[INDEX2(7,4,8)]+=tmp0_1;                              EM_S[INDEX2(0,5,8)]+=D_p[0]*w18;
5187                              EM_S[INDEX2(0,5,8)]+=tmp0_1;                              EM_S[INDEX2(1,5,8)]+=D_p[0]*w18;
5188                              EM_S[INDEX2(1,5,8)]+=tmp0_1;                              EM_S[INDEX2(2,5,8)]+=D_p[0]*w18;
5189                              EM_S[INDEX2(2,5,8)]+=tmp0_1;                              EM_S[INDEX2(3,5,8)]+=D_p[0]*w18;
5190                              EM_S[INDEX2(3,5,8)]+=tmp0_1;                              EM_S[INDEX2(4,5,8)]+=D_p[0]*w18;
5191                              EM_S[INDEX2(4,5,8)]+=tmp0_1;                              EM_S[INDEX2(5,5,8)]+=D_p[0]*w18;
5192                              EM_S[INDEX2(5,5,8)]+=tmp0_1;                              EM_S[INDEX2(6,5,8)]+=D_p[0]*w18;
5193                              EM_S[INDEX2(6,5,8)]+=tmp0_1;                              EM_S[INDEX2(7,5,8)]+=D_p[0]*w18;
5194                              EM_S[INDEX2(7,5,8)]+=tmp0_1;                              EM_S[INDEX2(0,6,8)]+=D_p[0]*w18;
5195                              EM_S[INDEX2(0,6,8)]+=tmp0_1;                              EM_S[INDEX2(1,6,8)]+=D_p[0]*w18;
5196                              EM_S[INDEX2(1,6,8)]+=tmp0_1;                              EM_S[INDEX2(2,6,8)]+=D_p[0]*w18;
5197                              EM_S[INDEX2(2,6,8)]+=tmp0_1;                              EM_S[INDEX2(3,6,8)]+=D_p[0]*w18;
5198                              EM_S[INDEX2(3,6,8)]+=tmp0_1;                              EM_S[INDEX2(4,6,8)]+=D_p[0]*w18;
5199                              EM_S[INDEX2(4,6,8)]+=tmp0_1;                              EM_S[INDEX2(5,6,8)]+=D_p[0]*w18;
5200                              EM_S[INDEX2(5,6,8)]+=tmp0_1;                              EM_S[INDEX2(6,6,8)]+=D_p[0]*w18;
5201                              EM_S[INDEX2(6,6,8)]+=tmp0_1;                              EM_S[INDEX2(7,6,8)]+=D_p[0]*w18;
5202                              EM_S[INDEX2(7,6,8)]+=tmp0_1;                              EM_S[INDEX2(0,7,8)]+=D_p[0]*w18;
5203                              EM_S[INDEX2(0,7,8)]+=tmp0_1;                              EM_S[INDEX2(1,7,8)]+=D_p[0]*w18;
5204                              EM_S[INDEX2(1,7,8)]+=tmp0_1;                              EM_S[INDEX2(2,7,8)]+=D_p[0]*w18;
5205                              EM_S[INDEX2(2,7,8)]+=tmp0_1;                              EM_S[INDEX2(3,7,8)]+=D_p[0]*w18;
5206                              EM_S[INDEX2(3,7,8)]+=tmp0_1;                              EM_S[INDEX2(4,7,8)]+=D_p[0]*w18;
5207                              EM_S[INDEX2(4,7,8)]+=tmp0_1;                              EM_S[INDEX2(5,7,8)]+=D_p[0]*w18;
5208                              EM_S[INDEX2(5,7,8)]+=tmp0_1;                              EM_S[INDEX2(6,7,8)]+=D_p[0]*w18;
5209                              EM_S[INDEX2(6,7,8)]+=tmp0_1;                              EM_S[INDEX2(7,7,8)]+=D_p[0]*w18;
                             EM_S[INDEX2(7,7,8)]+=tmp0_1;  
5210                          }                          }
5211                          ///////////////                          ///////////////
5212                          // process X //                          // process X //
# Line 5281  void Brick::assemblePDESingleReduced(Pas Line 5214  void Brick::assemblePDESingleReduced(Pas
5214                          if (!X.isEmpty()) {                          if (!X.isEmpty()) {
5215                              add_EM_F=true;                              add_EM_F=true;
5216                              const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                              const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
5217                              const double X_0 = X_p[0];                              const double wX0 = 8*X_p[0]*w12;
5218                              const double X_1 = X_p[1];                              const double wX1 = 8*X_p[1]*w13;
5219                              const double X_2 = X_p[2];                              const double wX2 = 8*X_p[2]*w14;
5220                              const double tmp0_1 = X_2*w21;                              EM_F[0]+=-wX0 - wX1 - wX2;
5221                              const double tmp1_1 = X_0*w19;                              EM_F[1]+= wX0 - wX1 - wX2;
5222                              const double tmp2_1 = X_1*w20;                              EM_F[2]+=-wX0 + wX1 - wX2;
5223                              const double tmp3_1 = X_0*w22;                              EM_F[3]+= wX0 + wX1 - wX2;
5224                              const double tmp4_1 = X_1*w23;                              EM_F[4]+=-wX0 - wX1 + wX2;
5225                              const double tmp5_1 = X_2*w24;                              EM_F[5]+= wX0 - wX1 + wX2;
5226                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_F[6]+=-wX0 + wX1 + wX2;
5227                              EM_F[1]+=tmp0_1 + tmp2_1 + tmp3_1;                              EM_F[7]+= wX0 + wX1 + wX2;
                             EM_F[2]+=tmp0_1 + tmp1_1 + tmp4_1;  
                             EM_F[3]+=tmp0_1 + tmp3_1 + tmp4_1;  
                             EM_F[4]+=tmp1_1 + tmp2_1 + tmp5_1;  
                             EM_F[5]+=tmp2_1 + tmp3_1 + tmp5_1;  
                             EM_F[6]+=tmp1_1 + tmp4_1 + tmp5_1;  
                             EM_F[7]+=tmp3_1 + tmp4_1 + tmp5_1;  
5228                          }                          }
5229                          ///////////////                          ///////////////
5230                          // process Y //                          // process Y //
# Line 5305  void Brick::assemblePDESingleReduced(Pas Line 5232  void Brick::assemblePDESingleReduced(Pas
5232                          if (!Y.isEmpty()) {                          if (!Y.isEmpty()) {
5233                              add_EM_F=true;                              add_EM_F=true;
5234                              const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                              const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
5235                              const double tmp0_1 = Y_p[0]*w25;                              EM_F[0]+=8*Y_p[0]*w18;
5236                              EM_F[0]+=tmp0_1;                              EM_F[1]+=8*Y_p[0]*w18;
5237                              EM_F[1]+=tmp0_1;                              EM_F[2]+=8*Y_p[0]*w18;
5238                              EM_F[2]+=tmp0_1;                              EM_F[3]+=8*Y_p[0]*w18;
5239                              EM_F[3]+=tmp0_1;                              EM_F[4]+=8*Y_p[0]*w18;
5240                              EM_F[4]+=tmp0_1;                              EM_F[5]+=8*Y_p[0]*w18;
5241                              EM_F[5]+=tmp0_1;                              EM_F[6]+=8*Y_p[0]*w18;
5242                              EM_F[6]+=tmp0_1;                              EM_F[7]+=8*Y_p[0]*w18;
                             EM_F[7]+=tmp0_1;  
5243                          }                          }
5244    
5245                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)
# Line 7472  void Brick::assemblePDESystemReduced(Pas Line 7398  void Brick::assemblePDESystemReduced(Pas
7398          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
7399          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
7400  {  {
     const double h0 = m_dx[0];  
     const double h1 = m_dx[1];  
     const double h2 = m_dx[2];  
7401      dim_t numEq, numComp;      dim_t numEq, numComp;
7402      if (!mat)      if (!mat)
7403          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
# Line 7483  void Brick::assemblePDESystemReduced(Pas Line 7406  void Brick::assemblePDESystemReduced(Pas
7406          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
7407      }      }
7408    
7409      const double w0 = .0625*h0;      const double w0 = .0625*m_dx[0];
7410      const double w1 = .0625*h1;      const double w1 = .0625*m_dx[1];
7411      const double w2 = .0625*h2;      const double w2 = .0625*m_dx[2];
7412      const double w3 = .03125*h0*h1;      const double w3 = .03125*m_dx[0]*m_dx[1];
7413      const double w4 = .03125*h0*h2;      const double w4 = .03125*m_dx[0]*m_dx[2];
7414      const double w5 = .03125*h1*h2;      const double w5 = .03125*m_dx[1]*m_dx[2];
7415      const double w6 = .0625*h0*h1/h2;      const double w6 = .0625*m_dx[0]*m_dx[1]/m_dx[2];
7416      const double w7 = .0625*h0*h2/h1;      const double w7 = .0625*m_dx[0]*m_dx[2]/m_dx[1];
7417      const double w8 = .0625*h1*h2/h0;      const double w8 = .0625*m_dx[1]*m_dx[2]/m_dx[0];
7418      const double w9 = .015625*h0*h1*h2;      const double w9 = .015625*m_dx[0]*m_dx[1]*m_dx[2];
7419    
7420      rhs.requireWrite();      rhs.requireWrite();
7421  #pragma omp parallel  #pragma omp parallel
# Line 7692  void Brick::assemblePDESystemReduced(Pas Line 7615  void Brick::assemblePDESystemReduced(Pas
7615                                      const double tmp0 = C_p[INDEX3(k, m, 0, numEq, numComp)]*w5;                                      const double tmp0 = C_p[INDEX3(k, m, 0, numEq, numComp)]*w5;
7616                                      const double tmp1 = C_p[INDEX3(k, m, 1, numEq, numComp)]*w4;                                      const double tmp1 = C_p[INDEX3(k, m, 1, numEq, numComp)]*w4;
7617                                      const double tmp2 = C_p[INDEX3(k, m, 2, numEq, numComp)]*w3;                                      const double tmp2 = C_p[INDEX3(k, m, 2, numEq, numComp)]*w3;
7618                                      const double tmp2_1 = tmp0;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;
7619                                      const double tmp5_1 = -tmp0;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;
7620                                      const double tmp0_1 = tmp1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;
7621                                      const double tmp4_1 = -tmp1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;
7622                                      const double tmp3_1 = tmp2;                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;
7623                                      const double tmp1_1 = -tmp2;                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;
7624                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;
7625                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;
7626                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;
7627                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;
7628                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;
7629                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;
7630                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;
7631                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;
7632                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;
7633                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;
7634                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;
7635                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;
7636                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;
7637                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;
7638                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;
7639                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;
7640                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;
7641                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;
7642                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;
7643                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;
7644                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;
7645                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;
7646                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;
7647                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;
7648                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;
7649                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;
7650                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;
7651                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;
7652                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;
7653                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;
7654                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;
7655                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;
7656                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;
7657                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;
7658                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;
7659                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;
7660                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;
7661                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;
7662                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;
7663                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;
7664                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;
7665                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;
7666                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;
7667                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;
7668                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;
7669                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;
7670                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;
7671                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;
7672                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;
7673                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;
7674                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;
7675                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;
7676                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;
7677                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;
7678                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;
7679                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;
7680                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;
7681                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;
                                     EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;  
7682                                  }                                  }
7683                              }                              }
7684                          }                          }
# Line 7773  void Brick::assemblePDESystemReduced(Pas Line 7690  void Brick::assemblePDESystemReduced(Pas
7690                              const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                              const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
7691                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
7692                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
7693                                      const double D_0 = D_p[INDEX2(k, m, numEq)];                                      const double tmp0 = D_p[INDEX2(k, m, numEq)]*w9;
7694                                      const double tmp0_1 = D_0*w9;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0;
7695                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0;
7696                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0;
7697                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp0;
7698                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0;
7699                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp0;
7700                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp0;
7701                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=tmp0;
7702                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0;
7703                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0;
7704                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0;
7705                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0;
7706                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp0;
7707                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0;
7708                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=tmp0;
7709                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp0;
7710                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0;
7711                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0;
7712                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0;
7713                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0;
7714                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp0;
7715                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=tmp0;
7716                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0;
7717                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0;
7718                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp0;
7719                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0;
7720                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0;
7721                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0;
7722                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=tmp0;
7723                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp0;
7724                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0;
7725                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0;
7726                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0;
7727                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp0;
7728                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp0;
7729                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=tmp0;
7730                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0;
7731                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0;
7732                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0;
7733                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0;
7734                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp0;
7735                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0;
7736                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=tmp0;
7737                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp0;
7738                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0;
7739                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0;
7740                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0;
7741                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0;
7742                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp0;
7743                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=tmp0;
7744                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0;
7745                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0;
7746                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0;
7747                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0;
7748                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0;
7749                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0;
7750                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=tmp0;
7751                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp0;
7752                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0;
7753                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0;
7754                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0;
7755                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0;
7756                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0;
7757                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=tmp0_1;                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0;
                                     EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp0_1;  
7758                                  }                                  }
7759                              }                              }
7760                          }                          }
# Line 7849  void Brick::assemblePDESystemReduced(Pas Line 7765  void Brick::assemblePDESystemReduced(Pas
7765                              add_EM_F=true;                              add_EM_F=true;
7766                              const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                              const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
7767                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
7768                                  const double X_0 = X_p[INDEX2(k, 0, numEq)];                                  const double tmp0 = 8*X_p[INDEX2(k, 0, numEq)]*w5;
7769                                  const double X_1 = X_p[INDEX2(k, 1, numEq)];                                  const double tmp1 = 8*X_p[INDEX2(k, 1, numEq)]*w4;
7770                                  const double X_2 = X_p[INDEX2(k, 2, numEq)];                                  const double tmp2 = 8*X_p[INDEX2(k, 2, numEq)]*w3;
7771                                  const double tmp1_1 = -8*X_0*w5;                                  EM_F[INDEX2(k,0,numEq)]+=-tmp0 - tmp1 - tmp2;
7772                                  const double tmp2_1 = -8*X_1*w4;                                  EM_F[INDEX2(k,1,numEq)]+= tmp0 - tmp1 - tmp2;
7773                                  const double tmp3_1 = 8*X_0*w5;                                  EM_F[INDEX2(k,2,numEq)]+=-tmp0 + tmp1 - tmp2;
7774                                  const double tmp4_1 = 8*X_1*w4;                                  EM_F[INDEX2(k,3,numEq)]+= tmp0 + tmp1 - tmp2;
7775                                  const double tmp5_1 = 8*X_2*w3;                                  EM_F[INDEX2(k,4,numEq)]+=-tmp0 - tmp1 + tmp2;
7776                                  const double tmp0_1 = -8*X_2*w3;                                  EM_F[INDEX2(k,5,numEq)]+= tmp0 - tmp1 + tmp2;
7777                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                  EM_F[INDEX2(k,6,numEq)]+=-tmp0 + tmp1 + tmp2;
7778                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1 + tmp3_1;                                  EM_F[INDEX2(k,7,numEq)]+= tmp0 + tmp1 + tmp2;
                                 EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp1_1 + tmp4_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp0_1 + tmp3_1 + tmp4_1;  
                                 EM_F[INDEX2(k,4,numEq)]+=tmp1_1 + tmp2_1 + tmp5_1;  
                                 EM_F[INDEX2(k,5,numEq)]+=tmp2_1 + tmp3_1 + tmp5_1;  
                                 EM_F[INDEX2(k,6,numEq)]+=tmp1_1 + tmp4_1 + tmp5_1;  
                                 EM_F[INDEX2(k,7,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;  
7779                              }                              }
7780                          }                          }
7781                          ///////////////                          ///////////////
# Line 7875  void Brick::assemblePDESystemReduced(Pas Line 7785  void Brick::assemblePDESystemReduced(Pas
7785                              add_EM_F=true;                              add_EM_F=true;
7786                              const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                              const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
7787                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
7788                                  const double Y_0 = Y_p[k];                                  EM_F[INDEX2(k,0,numEq)]+=8*Y_p[k]*w9;
7789                                  const double tmp0_1 = 8*Y_0*w9;                                  EM_F[INDEX2(k,1,numEq)]+=8*Y_p[k]*w9;
7790                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,2,numEq)]+=8*Y_p[k]*w9;
7791                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,3,numEq)]+=8*Y_p[k]*w9;
7792                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,4,numEq)]+=8*Y_p[k]*w9;
7793                                  EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,5,numEq)]+=8*Y_p[k]*w9;
7794                                  EM_F[INDEX2(k,4,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,6,numEq)]+=8*Y_p[k]*w9;
7795                                  EM_F[INDEX2(k,5,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,7,numEq)]+=8*Y_p[k]*w9;
                                 EM_F[INDEX2(k,6,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,7,numEq)]+=tmp0_1;  
7796                              }                              }
7797                          }                          }
7798    

Legend:
Removed from v.4370  
changed lines
  Added in v.4375

  ViewVC Help
Powered by ViewVC 1.1.26