/[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 4375 by caltinay, Mon Apr 22 05:35:52 2013 UTC revision 4377 by caltinay, Tue Apr 23 05:30:02 2013 UTC
# Line 7438  void Brick::assemblePDESystemReduced(Pas Line 7438  void Brick::assemblePDESystemReduced(Pas
7438                              const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                              const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
7439                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
7440                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
7441                                      const double A_00 = A_p[INDEX4(k,0,m,0,numEq,3,numComp)];                                      const double Aw00 = A_p[INDEX4(k,0,m,0,numEq,3,numComp)]*w8;
7442                                      const double A_01 = A_p[INDEX4(k,0,m,1,numEq,3,numComp)];                                      const double Aw10 = A_p[INDEX4(k,1,m,0,numEq,3,numComp)]*w2;
7443                                      const double A_02 = A_p[INDEX4(k,0,m,2,numEq,3,numComp)];                                      const double Aw20 = A_p[INDEX4(k,2,m,0,numEq,3,numComp)]*w1;
7444                                      const double A_10 = A_p[INDEX4(k,1,m,0,numEq,3,numComp)];                                      const double Aw01 = A_p[INDEX4(k,0,m,1,numEq,3,numComp)]*w2;
7445                                      const double A_11 = A_p[INDEX4(k,1,m,1,numEq,3,numComp)];                                      const double Aw11 = A_p[INDEX4(k,1,m,1,numEq,3,numComp)]*w7;
7446                                      const double A_12 = A_p[INDEX4(k,1,m,2,numEq,3,numComp)];                                      const double Aw21 = A_p[INDEX4(k,2,m,1,numEq,3,numComp)]*w0;
7447                                      const double A_20 = A_p[INDEX4(k,2,m,0,numEq,3,numComp)];                                      const double Aw02 = A_p[INDEX4(k,0,m,2,numEq,3,numComp)]*w1;
7448                                      const double A_21 = A_p[INDEX4(k,2,m,1,numEq,3,numComp)];                                      const double Aw12 = A_p[INDEX4(k,1,m,2,numEq,3,numComp)]*w0;
7449                                      const double A_22 = A_p[INDEX4(k,2,m,2,numEq,3,numComp)];                                      const double Aw22 = A_p[INDEX4(k,2,m,2,numEq,3,numComp)]*w6;
7450                                      const double tmp0 = (A_01 + A_10)*w2;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+= Aw00 + Aw01 + Aw02 + Aw10 + Aw11 + Aw12 + Aw20 + Aw21 + Aw22;
7451                                      const double tmp1 = (A_02 + A_20)*w1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=-Aw00 - Aw01 - Aw02 + Aw10 + Aw11 + Aw12 + Aw20 + Aw21 + Aw22;
7452                                      const double tmp2 = (A_12 + A_21)*w0;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+= Aw00 + Aw01 + Aw02 - Aw10 - Aw11 - Aw12 + Aw20 + Aw21 + Aw22;
7453                                      const double tmp3 = A_00*w8;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=-Aw00 - Aw01 - Aw02 - Aw10 - Aw11 - Aw12 + Aw20 + Aw21 + Aw22;
7454                                      const double tmp4 = A_01*w2;                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+= Aw00 + Aw01 + Aw02 + Aw10 + Aw11 + Aw12 - Aw20 - Aw21 - Aw22;
7455                                      const double tmp5 = A_02*w1;                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=-Aw00 - Aw01 - Aw02 + Aw10 + Aw11 + Aw12 - Aw20 - Aw21 - Aw22;
7456                                      const double tmp6 = A_10*w2;                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+= Aw00 + Aw01 + Aw02 - Aw10 - Aw11 - Aw12 - Aw20 - Aw21 - Aw22;
7457                                      const double tmp7 = A_11*w7;                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=-Aw00 - Aw01 - Aw02 - Aw10 - Aw11 - Aw12 - Aw20 - Aw21 - Aw22;
7458                                      const double tmp8 = A_12*w0;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=-Aw00 + Aw01 + Aw02 - Aw10 + Aw11 + Aw12 - Aw20 + Aw21 + Aw22;
7459                                      const double tmp9 = A_20*w1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+= Aw00 - Aw01 - Aw02 - Aw10 + Aw11 + Aw12 - Aw20 + Aw21 + Aw22;
7460                                      const double tmp10 = A_21*w0;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=-Aw00 + Aw01 + Aw02 + Aw10 - Aw11 - Aw12 - Aw20 + Aw21 + Aw22;
7461                                      const double tmp11 = A_22*w6;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+= Aw00 - Aw01 - Aw02 + Aw10 - Aw11 - Aw12 - Aw20 + Aw21 + Aw22;
7462                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp7 + tmp11;                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=-Aw00 + Aw01 + Aw02 - Aw10 + Aw11 + Aw12 + Aw20 - Aw21 - Aw22;
7463                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp2 - tmp3 - tmp4 - tmp5 + tmp6 + tmp7 + tmp9 + tmp11;                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+= Aw00 - Aw01 - Aw02 - Aw10 + Aw11 + Aw12 + Aw20 - Aw21 - Aw22;
7464                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp1 + tmp3 + tmp4 - tmp6 - tmp7 - tmp8 + tmp10 + tmp11;                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=-Aw00 + Aw01 + Aw02 + Aw10 - Aw11 - Aw12 + Aw20 - Aw21 - Aw22;
7465                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=-tmp0 - tmp7 - tmp8 + tmp10 - tmp5 + tmp9 - tmp3 + tmp11;                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+= Aw00 - Aw01 - Aw02 + Aw10 - Aw11 - Aw12 + Aw20 - Aw21 - Aw22;
7466                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0 + tmp3 + tmp5 + tmp7 + tmp8 - tmp9 - tmp10 - tmp11;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+= Aw00 - Aw01 + Aw02 + Aw10 - Aw11 + Aw12 + Aw20 - Aw21 + Aw22;
7467                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=-tmp10 + tmp8 - tmp1 + tmp6 - tmp4 - tmp11 + tmp7 - tmp3;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=-Aw00 + Aw01 - Aw02 + Aw10 - Aw11 + Aw12 + Aw20 - Aw21 + Aw22;
7468                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=-tmp7 + tmp4 - tmp6 - tmp2 - tmp9 + tmp5 + tmp3 - tmp11;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+= Aw00 - Aw01 + Aw02 - Aw10 + Aw11 - Aw12 + Aw20 - Aw21 + Aw22;
7469                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp7 - tmp2 - tmp11 - tmp3;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=-Aw00 + Aw01 - Aw02 - Aw10 + Aw11 - Aw12 + Aw20 - Aw21 + Aw22;
7470                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp4 - tmp6 - tmp9 + tmp5 + tmp2 + tmp7 - tmp3 + tmp11;                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+= Aw00 - Aw01 + Aw02 + Aw10 - Aw11 + Aw12 - Aw20 + Aw21 - Aw22;
7471                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2 + tmp3 + tmp7 + tmp11;                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=-Aw00 + Aw01 - Aw02 + Aw10 - Aw11 + Aw12 - Aw20 + Aw21 - Aw22;
7472                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0 - tmp7 - tmp9 - tmp8 + tmp5 + tmp10 - tmp3 + tmp11;                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+= Aw00 - Aw01 + Aw02 - Aw10 + Aw11 - Aw12 - Aw20 + Aw21 - Aw22;
7473                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=-tmp7 - tmp1 + tmp6 - tmp8 - tmp4 + tmp3 + tmp10 + tmp11;                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=-Aw00 + Aw01 - Aw02 - Aw10 + Aw11 - Aw12 - Aw20 + Aw21 - Aw22;
7474                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp4 - tmp10 + tmp8 - tmp6 + tmp1 - tmp11 + tmp7 - tmp3;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=-Aw00 - Aw01 + Aw02 - Aw10 - Aw11 + Aw12 - Aw20 - Aw21 + Aw22;
7475                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=-tmp0 - tmp10 + tmp8 + tmp3 - tmp11 + tmp7 - tmp5 + tmp9;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+= Aw00 + Aw01 - Aw02 - Aw10 - Aw11 + Aw12 - Aw20 - Aw21 + Aw22;
7476                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=tmp0 - tmp7 - tmp2 + tmp1 - tmp11 - tmp3;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=-Aw00 - Aw01 + Aw02 + Aw10 + Aw11 - Aw12 - Aw20 - Aw21 + Aw22;
7477                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=-tmp7 + tmp6 - tmp2 - tmp4 + tmp3 - tmp11 - tmp5 + tmp9;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+= Aw00 + Aw01 - Aw02 + Aw10 + Aw11 - Aw12 - Aw20 - Aw21 + Aw22;
7478                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=-tmp7 - tmp10 + tmp8 + tmp6 - tmp4 + tmp1 + tmp3 + tmp11;                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=-Aw00 - Aw01 + Aw02 - Aw10 - Aw11 + Aw12 + Aw20 + Aw21 - Aw22;
7479                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0 - tmp7 - tmp10 + tmp8 - tmp5 + tmp9 - tmp3 + tmp11;                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+= Aw00 + Aw01 - Aw02 - Aw10 - Aw11 + Aw12 + Aw20 + Aw21 - Aw22;
7480                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2 + tmp3 + tmp7 + tmp11;                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=-Aw00 - Aw01 + Aw02 + Aw10 + Aw11 - Aw12 + Aw20 + Aw21 - Aw22;
7481                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp4 - tmp6 - tmp2 + tmp7 - tmp5 + tmp9 - tmp3 + tmp11;                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+= Aw00 + Aw01 - Aw02 + Aw10 + Aw11 - Aw12 + Aw20 + Aw21 - Aw22;
7482                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=-tmp7 + tmp6 - tmp9 - tmp4 + tmp5 + tmp2 + tmp3 - tmp11;                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+= Aw00 + Aw01 - Aw02 + Aw10 + Aw11 - Aw12 + Aw20 + Aw21 - Aw22;
7483                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=tmp0 - tmp1 + tmp2 - tmp3 - tmp7 - tmp11;                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=-Aw00 - Aw01 + Aw02 + Aw10 + Aw11 - Aw12 + Aw20 + Aw21 - Aw22;
7484                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=-tmp0 + tmp3 + tmp5 + tmp7 - tmp8 - tmp9 + tmp10 - tmp11;                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+= Aw00 + Aw01 - Aw02 - Aw10 - Aw11 + Aw12 + Aw20 + Aw21 - Aw22;
7485                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=-tmp1 - tmp3 + tmp4 - tmp6 + tmp7 - tmp8 + tmp10 - tmp11;                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=-Aw00 - Aw01 + Aw02 - Aw10 - Aw11 + Aw12 + Aw20 + Aw21 - Aw22;
7486                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=-tmp0 - tmp7 + tmp8 - tmp9 + tmp5 - tmp3 - tmp10 + tmp11;                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+= Aw00 + Aw01 - Aw02 + Aw10 + Aw11 - Aw12 - Aw20 - Aw21 + Aw22;
7487                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=-tmp1 + tmp3 + tmp4 - tmp6 - tmp7 + tmp8 - tmp10 + tmp11;                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=-Aw00 - Aw01 + Aw02 + Aw10 + Aw11 - Aw12 - Aw20 - Aw21 + Aw22;
7488                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp6 - tmp2 - tmp3 - tmp4 + tmp5 + tmp7 - tmp9 + tmp11;                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+= Aw00 + Aw01 - Aw02 - Aw10 - Aw11 + Aw12 - Aw20 - Aw21 + Aw22;
7489                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0 - tmp1 - tmp2 + tmp3 + tmp7 + tmp11;                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=-Aw00 - Aw01 + Aw02 - Aw10 - Aw11 + Aw12 - Aw20 - Aw21 + Aw22;
7490                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2 - tmp3 - tmp7 - tmp11;                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=-Aw00 + Aw01 - Aw02 - Aw10 + Aw11 - Aw12 - Aw20 + Aw21 - Aw22;
7491                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=-tmp7 + tmp4 - tmp6 + tmp2 + tmp3 - tmp11 - tmp5 + tmp9;                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+= Aw00 - Aw01 + Aw02 - Aw10 + Aw11 - Aw12 - Aw20 + Aw21 - Aw22;
7492                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp1 - tmp3 - tmp4 + tmp6 + tmp7 - tmp8 + tmp10 - tmp11;                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=-Aw00 + Aw01 - Aw02 + Aw10 - Aw11 + Aw12 - Aw20 + Aw21 - Aw22;
7493                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0 - tmp8 + tmp3 - tmp11 + tmp7 + tmp10 - tmp5 + tmp9;                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+= Aw00 - Aw01 + Aw02 + Aw10 - Aw11 + Aw12 - Aw20 + Aw21 - Aw22;
7494                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0 - tmp8 + tmp3 - tmp11 + tmp7 + tmp10 - tmp5 + tmp9;                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=-Aw00 + Aw01 - Aw02 - Aw10 + Aw11 - Aw12 + Aw20 - Aw21 + Aw22;
7495                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp6 - tmp8 - tmp4 + tmp1 - tmp11 + tmp7 + tmp10 - tmp3;                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+= Aw00 - Aw01 + Aw02 - Aw10 + Aw11 - Aw12 + Aw20 - Aw21 + Aw22;
7496                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=-tmp7 + tmp4 - tmp6 + tmp2 + tmp3 - tmp11 - tmp5 + tmp9;                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=-Aw00 + Aw01 - Aw02 + Aw10 - Aw11 + Aw12 + Aw20 - Aw21 + Aw22;
7497                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=-tmp0 - tmp7 + tmp1 + tmp2 - tmp11 - tmp3;                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+= Aw00 - Aw01 + Aw02 + Aw10 - Aw11 + Aw12 + Aw20 - Aw21 + Aw22;
7498                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0 - tmp1 - tmp2 + tmp3 + tmp7 + tmp11;                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+= Aw00 - Aw01 - Aw02 + Aw10 - Aw11 - Aw12 + Aw20 - Aw21 - Aw22;
7499                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp6 - tmp2 - tmp9 - tmp4 + tmp5 + tmp7 - tmp3 + tmp11;                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=-Aw00 + Aw01 + Aw02 + Aw10 - Aw11 - Aw12 + Aw20 - Aw21 - Aw22;
7500                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=-tmp7 + tmp4 - tmp10 + tmp8 - tmp1 - tmp6 + tmp3 + tmp11;                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+= Aw00 - Aw01 - Aw02 - Aw10 + Aw11 + Aw12 + Aw20 - Aw21 - Aw22;
7501                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=-tmp0 - tmp7 - tmp10 + tmp8 - tmp9 + tmp5 - tmp3 + tmp11;                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=-Aw00 + Aw01 + Aw02 - Aw10 + Aw11 + Aw12 + Aw20 - Aw21 - Aw22;
7502                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp4 - tmp1 - tmp6 - tmp8 - tmp11 + tmp7 + tmp10 - tmp3;                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+= Aw00 - Aw01 - Aw02 + Aw10 - Aw11 - Aw12 - Aw20 + Aw21 + Aw22;
7503                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=-tmp0 - tmp9 - tmp8 + tmp5 + tmp3 - tmp11 + tmp7 + tmp10;                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=-Aw00 + Aw01 + Aw02 + Aw10 - Aw11 - Aw12 - Aw20 + Aw21 + Aw22;
7504                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=tmp0 - tmp7 - tmp1 + tmp2 - tmp11 - tmp3;                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+= Aw00 - Aw01 - Aw02 - Aw10 + Aw11 + Aw12 - Aw20 + Aw21 + Aw22;
7505                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=-tmp7 + tmp6 - tmp9 - tmp4 + tmp5 + tmp2 + tmp3 - tmp11;                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=-Aw00 + Aw01 + Aw02 - Aw10 + Aw11 + Aw12 - Aw20 + Aw21 + Aw22;
7506                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp4 - tmp6 - tmp2 + tmp7 - tmp5 + tmp9 - tmp3 + tmp11;                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=-Aw00 - Aw01 - Aw02 - Aw10 - Aw11 - Aw12 - Aw20 - Aw21 - Aw22;
7507                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=-tmp0 - tmp2 + tmp1 + tmp3 + tmp7 + tmp11;                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+= Aw00 + Aw01 + Aw02 - Aw10 - Aw11 - Aw12 - Aw20 - Aw21 - Aw22;
7508                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0 - tmp7 - tmp10 + tmp8 - tmp5 + tmp9 - tmp3 + tmp11;                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=-Aw00 - Aw01 - Aw02 + Aw10 + Aw11 + Aw12 - Aw20 - Aw21 - Aw22;
7509                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=-tmp7 - tmp10 + tmp8 + tmp6 - tmp4 + tmp1 + tmp3 + tmp11;                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+= Aw00 + Aw01 + Aw02 + Aw10 + Aw11 + Aw12 - Aw20 - Aw21 - Aw22;
7510                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=-tmp7 + tmp6 - tmp2 - tmp4 + tmp3 - tmp11 - tmp5 + tmp9;                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=-Aw00 - Aw01 - Aw02 - Aw10 - Aw11 - Aw12 + Aw20 + Aw21 + Aw22;
7511                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=tmp0 - tmp7 - tmp2 + tmp1 - tmp11 - tmp3;                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+= Aw00 + Aw01 + Aw02 - Aw10 - Aw11 - Aw12 + Aw20 + Aw21 + Aw22;
7512                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=-tmp0 - tmp10 + tmp8 + tmp3 - tmp11 + tmp7 - tmp5 + tmp9;                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=-Aw00 - Aw01 - Aw02 + Aw10 + Aw11 + Aw12 + Aw20 + Aw21 + Aw22;
7513                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp4 - tmp10 + tmp8 - tmp6 + tmp1 - tmp11 + tmp7 - tmp3;                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+= Aw00 + Aw01 + Aw02 + Aw10 + Aw11 + Aw12 + Aw20 + Aw21 + Aw22;
                                     EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=-tmp7 - tmp1 + tmp6 - tmp8 - tmp4 + tmp3 + tmp10 + tmp11;  
                                     EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0 - tmp7 - tmp9 - tmp8 + tmp5 + tmp10 - tmp3 + tmp11;  
                                     EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2 + tmp3 + tmp7 + tmp11;  
                                     EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp4 - tmp6 - tmp9 + tmp5 + tmp2 + tmp7 - tmp3 + tmp11;  
                                     EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=-tmp0 - tmp7 - tmp1 - tmp2 - tmp11 - tmp3;  
                                     EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=-tmp7 + tmp4 - tmp6 - tmp2 - tmp9 + tmp5 + tmp3 - tmp11;  
                                     EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=-tmp10 + tmp8 - tmp1 + tmp6 - tmp4 - tmp11 + tmp7 - tmp3;  
                                     EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0 + tmp3 + tmp5 + tmp7 + tmp8 - tmp9 - tmp10 - tmp11;  
                                     EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=-tmp0 - tmp3 - tmp5 - tmp7 - tmp8 + tmp9 + tmp10 + tmp11;  
                                     EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp1 + tmp3 + tmp4 - tmp6 - tmp7 - tmp8 + tmp10 + tmp11;  
                                     EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp2 - tmp3 - tmp4 - tmp5 + tmp6 + tmp7 + tmp9 + tmp11;  
                                     EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp7 + tmp11;  
7514                                  }                                  }
7515                              }                              }
7516                          }                          }
# Line 7534  void Brick::assemblePDESystemReduced(Pas Line 7522  void Brick::assemblePDESystemReduced(Pas
7522                              const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                              const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
7523                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
7524                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
7525                                      const double tmp0 = B_p[INDEX3(k,0,m, numEq, 3)]*w5;                                      const double wB0 = B_p[INDEX3(k,0,m, numEq, 3)]*w5;
7526                                      const double tmp1 = B_p[INDEX3(k,1,m, numEq, 3)]*w4;                                      const double wB1 = B_p[INDEX3(k,1,m, numEq, 3)]*w4;
7527                                      const double tmp2 = B_p[INDEX3(k,2,m, numEq, 3)]*w3;                                      const double wB2 = B_p[INDEX3(k,2,m, numEq, 3)]*w3;
7528                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=-wB0 - wB1 - wB2;
7529                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=-wB0 - wB1 - wB2;
7530                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=-wB0 - wB1 - wB2;
7531                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=-wB0 - wB1 - wB2;
7532                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=-wB0 - wB1 - wB2;
7533                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=-wB0 - wB1 - wB2;
7534                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=-wB0 - wB1 - wB2;
7535                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=-wB0 - wB1 - wB2;
7536                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+= wB0 - wB1 - wB2;
7537                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+= wB0 - wB1 - wB2;
7538                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+= wB0 - wB1 - wB2;
7539                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+= wB0 - wB1 - wB2;
7540                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+= wB0 - wB1 - wB2;
7541                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+= wB0 - wB1 - wB2;
7542                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+= wB0 - wB1 - wB2;
7543                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+= wB0 - wB1 - wB2;
7544                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=-wB0 + wB1 - wB2;
7545                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=-wB0 + wB1 - wB2;
7546                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=-wB0 + wB1 - wB2;
7547                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=-wB0 + wB1 - wB2;
7548                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=-wB0 + wB1 - wB2;
7549                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=-wB0 + wB1 - wB2;
7550                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=-wB0 + wB1 - wB2;
7551                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=-wB0 + wB1 - wB2;
7552                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+= wB0 + wB1 - wB2;
7553                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+= wB0 + wB1 - wB2;
7554                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+= wB0 + wB1 - wB2;
7555                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+= wB0 + wB1 - wB2;
7556                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+= wB0 + wB1 - wB2;
7557                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+= wB0 + wB1 - wB2;
7558                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+= wB0 + wB1 - wB2;
7559                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+= wB0 + wB1 - wB2;
7560                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=-wB0 - wB1 + wB2;
7561                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=-wB0 - wB1 + wB2;
7562                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=-wB0 - wB1 + wB2;
7563                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=-wB0 - wB1 + wB2;
7564                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=-wB0 - wB1 + wB2;
7565                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=-wB0 - wB1 + wB2;
7566                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=-wB0 - wB1 + wB2;
7567                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=-wB0 - wB1 + wB2;
7568                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+= wB0 - wB1 + wB2;
7569                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+= wB0 - wB1 + wB2;
7570                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+= wB0 - wB1 + wB2;
7571                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+= wB0 - wB1 + wB2;
7572                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+= wB0 - wB1 + wB2;
7573                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+= wB0 - wB1 + wB2;
7574                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+= wB0 - wB1 + wB2;
7575                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+= wB0 - wB1 + wB2;
7576                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=-wB0 + wB1 + wB2;
7577                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=-wB0 + wB1 + wB2;
7578                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=-wB0 + wB1 + wB2;
7579                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=-wB0 + wB1 + wB2;
7580                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=-wB0 + wB1 + wB2;
7581                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=-wB0 + wB1 + wB2;
7582                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=-wB0 + wB1 + wB2;
7583                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=-wB0 + wB1 + wB2;
7584                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=-tmp1 - tmp0 - tmp2;                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+= wB0 + wB1 + wB2;
7585                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+= wB0 + wB1 + wB2;
7586                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+= wB0 + wB1 + wB2;
7587                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+= wB0 + wB1 + wB2;
7588                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+= wB0 + wB1 + wB2;
7589                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+= wB0 + wB1 + wB2;
7590                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+= wB0 + wB1 + wB2;
7591                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+= wB0 + wB1 + wB2;
7592                                  }                                  }
7593                              }                              }
7594                          }                          }
# Line 7612  void Brick::assemblePDESystemReduced(Pas Line 7600  void Brick::assemblePDESystemReduced(Pas
7600                              const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                              const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
7601                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
7602                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
7603                                      const double tmp0 = C_p[INDEX3(k, m, 0, numEq, numComp)]*w5;                                      const double wC0 = C_p[INDEX3(k, m, 0, numEq, numComp)]*w5;
7604                                      const double tmp1 = C_p[INDEX3(k, m, 1, numEq, numComp)]*w4;                                      const double wC1 = C_p[INDEX3(k, m, 1, numEq, numComp)]*w4;
7605                                      const double tmp2 = C_p[INDEX3(k, m, 2, numEq, numComp)]*w3;                                      const double wC2 = C_p[INDEX3(k, m, 2, numEq, numComp)]*w3;
7606                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=-wC0 - wC1 - wC2;
7607                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=-wC0 - wC1 - wC2;
7608                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=-wC0 - wC1 - wC2;
7609                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=-wC0 - wC1 - wC2;
7610                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=-wC0 - wC1 - wC2;
7611                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=-wC0 - wC1 - wC2;
7612                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=-wC0 - wC1 - wC2;
7613                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=-tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=-wC0 - wC1 - wC2;
7614                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+= wC0 - wC1 - wC2;
7615                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+= wC0 - wC1 - wC2;
7616                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+= wC0 - wC1 - wC2;
7617                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+= wC0 - wC1 - wC2;
7618                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+= wC0 - wC1 - wC2;
7619                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+= wC0 - wC1 - wC2;
7620                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+= wC0 - wC1 - wC2;
7621                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+= tmp0 - tmp1 - tmp2;                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+= wC0 - wC1 - wC2;
7622                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=-wC0 + wC1 - wC2;
7623                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=-wC0 + wC1 - wC2;
7624                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=-wC0 + wC1 - wC2;
7625                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=-wC0 + wC1 - wC2;
7626                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=-wC0 + wC1 - wC2;
7627                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=-wC0 + wC1 - wC2;
7628                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=-wC0 + wC1 - wC2;
7629                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=-tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=-wC0 + wC1 - wC2;
7630                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+= wC0 + wC1 - wC2;
7631                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+= wC0 + wC1 - wC2;
7632                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+= wC0 + wC1 - wC2;
7633                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+= wC0 + wC1 - wC2;
7634                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+= wC0 + wC1 - wC2;
7635                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+= wC0 + wC1 - wC2;
7636                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+= wC0 + wC1 - wC2;
7637                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+= tmp0 + tmp1 - tmp2;                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+= wC0 + wC1 - wC2;
7638                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=-wC0 - wC1 + wC2;
7639                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=-wC0 - wC1 + wC2;
7640                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=-wC0 - wC1 + wC2;
7641                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=-wC0 - wC1 + wC2;
7642                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=-wC0 - wC1 + wC2;
7643                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=-wC0 - wC1 + wC2;
7644                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=-wC0 - wC1 + wC2;
7645                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=-tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=-wC0 - wC1 + wC2;
7646                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+= wC0 - wC1 + wC2;
7647                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+= wC0 - wC1 + wC2;
7648                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+= wC0 - wC1 + wC2;
7649                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+= wC0 - wC1 + wC2;
7650                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+= wC0 - wC1 + wC2;
7651                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+= wC0 - wC1 + wC2;
7652                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+= wC0 - wC1 + wC2;
7653                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+= tmp0 - tmp1 + tmp2;                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+= wC0 - wC1 + wC2;
7654                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=-wC0 + wC1 + wC2;
7655                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=-wC0 + wC1 + wC2;
7656                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=-wC0 + wC1 + wC2;
7657                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=-wC0 + wC1 + wC2;
7658                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=-wC0 + wC1 + wC2;
7659                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=-wC0 + wC1 + wC2;
7660                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=-wC0 + wC1 + wC2;
7661                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=-tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=-wC0 + wC1 + wC2;
7662                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+= wC0 + wC1 + wC2;
7663                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+= wC0 + wC1 + wC2;
7664                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+= wC0 + wC1 + wC2;
7665                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+= wC0 + wC1 + wC2;
7666                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+= wC0 + wC1 + wC2;
7667                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+= wC0 + wC1 + wC2;
7668                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+= wC0 + wC1 + wC2;
7669                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+= tmp0 + tmp1 + tmp2;                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+= wC0 + wC1 + wC2;
7670                                  }                                  }
7671                              }                              }
7672                          }                          }
# Line 7690  void Brick::assemblePDESystemReduced(Pas Line 7678  void Brick::assemblePDESystemReduced(Pas
7678                              const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                              const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
7679                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
7680                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
7681                                      const double tmp0 = D_p[INDEX2(k, m, numEq)]*w9;                                      const double wD = D_p[INDEX2(k, m, numEq)]*w9;
7682                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=wD;
7683                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=wD;
7684                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=wD;
7685                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=wD;
7686                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=wD;
7687                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=wD;
7688                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=wD;
7689                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=wD;
7690                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=wD;
7691                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=wD;
7692                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=wD;
7693                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=wD;
7694                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=wD;
7695                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=wD;
7696                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=wD;
7697                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=wD;
7698                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=wD;
7699                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=wD;
7700                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=wD;
7701                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=wD;
7702                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=wD;
7703                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=wD;
7704                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=wD;
7705                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=wD;
7706                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=wD;
7707                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=wD;
7708                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=wD;
7709                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=wD;
7710                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=wD;
7711                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=wD;
7712                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=wD;
7713                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=wD;
7714                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=wD;
7715                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=wD;
7716                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=wD;
7717                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=wD;
7718                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=wD;
7719                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=wD;
7720                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=wD;
7721                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=wD;
7722                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=wD;
7723                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=wD;
7724                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=wD;
7725                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=wD;
7726                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=wD;
7727                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=wD;
7728                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=wD;
7729                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=wD;
7730                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=wD;
7731                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=wD;
7732                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=wD;
7733                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=wD;
7734                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=wD;
7735                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=wD;
7736                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=wD;
7737                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=wD;
7738                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=wD;
7739                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=wD;
7740                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=wD;
7741                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=wD;
7742                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=wD;
7743                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=wD;
7744                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=wD;
7745                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=wD;
7746                                  }                                  }
7747                              }                              }
7748                          }                          }
# Line 7765  void Brick::assemblePDESystemReduced(Pas Line 7753  void Brick::assemblePDESystemReduced(Pas
7753                              add_EM_F=true;                              add_EM_F=true;
7754                              const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                              const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
7755                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
7756                                  const double tmp0 = 8*X_p[INDEX2(k, 0, numEq)]*w5;                                  const double wX0 = 8*X_p[INDEX2(k, 0, numEq)]*w5;
7757                                  const double tmp1 = 8*X_p[INDEX2(k, 1, numEq)]*w4;                                  const double wX1 = 8*X_p[INDEX2(k, 1, numEq)]*w4;
7758                                  const double tmp2 = 8*X_p[INDEX2(k, 2, numEq)]*w3;                                  const double wX2 = 8*X_p[INDEX2(k, 2, numEq)]*w3;
7759                                  EM_F[INDEX2(k,0,numEq)]+=-tmp0 - tmp1 - tmp2;                                  EM_F[INDEX2(k,0,numEq)]+=-wX0 - wX1 - wX2;
7760                                  EM_F[INDEX2(k,1,numEq)]+= tmp0 - tmp1 - tmp2;                                  EM_F[INDEX2(k,1,numEq)]+= wX0 - wX1 - wX2;
7761                                  EM_F[INDEX2(k,2,numEq)]+=-tmp0 + tmp1 - tmp2;                                  EM_F[INDEX2(k,2,numEq)]+=-wX0 + wX1 - wX2;
7762                                  EM_F[INDEX2(k,3,numEq)]+= tmp0 + tmp1 - tmp2;                                  EM_F[INDEX2(k,3,numEq)]+= wX0 + wX1 - wX2;
7763                                  EM_F[INDEX2(k,4,numEq)]+=-tmp0 - tmp1 + tmp2;                                  EM_F[INDEX2(k,4,numEq)]+=-wX0 - wX1 + wX2;
7764                                  EM_F[INDEX2(k,5,numEq)]+= tmp0 - tmp1 + tmp2;                                  EM_F[INDEX2(k,5,numEq)]+= wX0 - wX1 + wX2;
7765                                  EM_F[INDEX2(k,6,numEq)]+=-tmp0 + tmp1 + tmp2;                                  EM_F[INDEX2(k,6,numEq)]+=-wX0 + wX1 + wX2;
7766                                  EM_F[INDEX2(k,7,numEq)]+= tmp0 + tmp1 + tmp2;                                  EM_F[INDEX2(k,7,numEq)]+= wX0 + wX1 + wX2;
7767                              }                              }
7768                          }                          }
7769                          ///////////////                          ///////////////
# Line 7811  void Brick::assemblePDESystemReduced(Pas Line 7799  void Brick::assemblePDESystemReduced(Pas
7799  void Brick::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,  void Brick::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
7800        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
7801  {  {
7802      const double h0 = m_dx[0];      const double SQRT3 = 1.73205080756887719318;
7803      const double h1 = m_dx[1];      const double w12 = m_dx[0]*m_dx[1]/144;
7804      const double h2 = m_dx[2];      const double w11 = w12*(SQRT3 + 2);
7805      const double w0 = 0.0018607582807716854616*h1*h2;      const double w14 = w12*(4*SQRT3 + 7);
7806      const double w1 = 0.025917019497006092316*h1*h2;      const double w13 = w12*(-4*SQRT3 + 7);
7807      const double w2 = 0.0069444444444444444444*h1*h2;      const double w10 = w12*(-SQRT3 + 2);
7808      const double w3 = 0.00049858867864229740201*h1*h2;      const double w7 = m_dx[0]*m_dx[2]/144;
7809      const double w4 = 0.09672363354357992482*h1*h2;      const double w6 = w7*(SQRT3 + 2);
7810      const double w5 = 0.055555555555555555556*h1*h2;      const double w8 = w7*(-4*SQRT3 + 7);
7811      const double w6 = 0.11111111111111111111*h1*h2;      const double w9 = w7*(4*SQRT3 + 7);
7812      const double w7 = 0.027777777777777777778*h1*h2;      const double w5 = w7*(-SQRT3 + 2);
7813      const double w8 = 0.1555021169820365539*h1*h2;      const double w2 = m_dx[1]*m_dx[2]/144;
7814      const double w9 = 0.041666666666666666667*h1*h2;      const double w1 = w2*(SQRT3 + 2);
7815      const double w10 = 0.01116454968463011277*h1*h2;      const double w4 = w2*(4*SQRT3 + 7);
7816      const double w11 = 0.25*h1*h2;      const double w3 = w2*(-4*SQRT3 + 7);
7817      const double w12 = 0.0018607582807716854616*h0*h2;      const double w0 = w2*(-SQRT3 + 2);
     const double w13 = 0.025917019497006092316*h0*h2;  
     const double w14 = 0.0069444444444444444444*h0*h2;  
     const double w15 = 0.00049858867864229740201*h0*h2;  
     const double w16 = 0.09672363354357992482*h0*h2;  
     const double w17 = 0.055555555555555555556*h0*h2;  
     const double w18 = 0.11111111111111111111*h0*h2;  
     const double w19 = 0.027777777777777777778*h0*h2;  
     const double w20 = 0.1555021169820365539*h0*h2;  
     const double w21 = 0.041666666666666666667*h0*h2;  
     const double w22 = 0.01116454968463011277*h0*h2;  
     const double w23 = 0.25*h0*h2;  
     const double w24 = 0.0018607582807716854616*h0*h1;  
     const double w25 = 0.025917019497006092316*h0*h1;  
     const double w26 = 0.0069444444444444444444*h0*h1;  
     const double w27 = 0.00049858867864229740201*h0*h1;  
     const double w28 = 0.09672363354357992482*h0*h1;  
     const double w29 = 0.055555555555555555556*h0*h1;  
     const double w30 = 0.027777777777777777778*h0*h1;  
     const double w31 = 0.11111111111111111111*h0*h1;  
     const double w32 = 0.1555021169820365539*h0*h1;  
     const double w33 = 0.041666666666666666667*h0*h1;  
     const double w34 = 0.01116454968463011277*h0*h1;  
     const double w35 = 0.25*h0*h1;  
7818      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
7819      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
7820      rhs.requireWrite();      rhs.requireWrite();
# Line 7873  void Brick::assemblePDEBoundarySingle(Pa Line 7838  void Brick::assemblePDEBoundarySingle(Pa
7838                                  const double d_1 = d_p[1];                                  const double d_1 = d_p[1];
7839                                  const double d_2 = d_p[2];                                  const double d_2 = d_p[2];
7840                                  const double d_3 = d_p[3];                                  const double d_3 = d_p[3];
7841                                  const double tmp3_0 = d_1 + d_3;                                  const double tmp0 = w0*(d_0 + d_1);
7842                                  const double tmp6_0 = d_0 + d_1 + d_2 + d_3;                                  const double tmp1 = w1*(d_2 + d_3);
7843                                  const double tmp2_0 = d_0 + d_2;                                  const double tmp2 = w0*(d_0 + d_2);
7844                                  const double tmp0_0 = d_0 + d_1;                                  const double tmp3 = w1*(d_1 + d_3);
7845                                  const double tmp4_0 = d_0 + d_3;                                  const double tmp4 = w0*(d_1 + d_3);
7846                                  const double tmp1_0 = d_2 + d_3;                                  const double tmp5 = w1*(d_0 + d_2);
7847                                  const double tmp5_0 = d_1 + d_2;                                  const double tmp6 = w0*(d_2 + d_3);
7848                                  const double tmp14_1 = d_0*w4;                                  const double tmp7 = w1*(d_0 + d_1);
7849                                  const double tmp17_1 = d_3*w4;                                  const double tmp8 = w2*(d_0 + d_3);
7850                                  const double tmp12_1 = d_1*w4;                                  const double tmp9 = w2*(d_1 + d_2);
7851                                  const double tmp9_1 = d_2*w4;                                  const double tmp10 = w2*(d_0 + d_1 + d_2 + d_3);
7852                                  const double tmp4_1 = tmp3_0*w0;                                  EM_S[INDEX2(0,0,8)]+=d_0*w4 + d_3*w3 + tmp9;
7853                                  const double tmp3_1 = tmp3_0*w1;                                  EM_S[INDEX2(0,2,8)]+=tmp6 + tmp7;
7854                                  const double tmp7_1 = tmp0_0*w1;                                  EM_S[INDEX2(0,4,8)]+=tmp4 + tmp5;
7855                                  const double tmp8_1 = d_1*w3;                                  EM_S[INDEX2(0,6,8)]+=tmp10;
7856                                  const double tmp16_1 = d_0*w3;                                  EM_S[INDEX2(2,0,8)]+=tmp6 + tmp7;
7857                                  const double tmp18_1 = tmp6_0*w2;                                  EM_S[INDEX2(2,2,8)]+=d_1*w4 + d_2*w3 + tmp8;
7858                                  const double tmp1_1 = tmp1_0*w1;                                  EM_S[INDEX2(2,4,8)]+=tmp10;
7859                                  const double tmp15_1 = tmp5_0*w2;                                  EM_S[INDEX2(2,6,8)]+=tmp2 + tmp3;
7860                                  const double tmp6_1 = tmp1_0*w0;                                  EM_S[INDEX2(4,0,8)]+=tmp4 + tmp5;
7861                                  const double tmp2_1 = tmp2_0*w0;                                  EM_S[INDEX2(4,2,8)]+=tmp10;
7862                                  const double tmp13_1 = d_3*w3;                                  EM_S[INDEX2(4,4,8)]+=d_1*w3 + d_2*w4 + tmp8;
7863                                  const double tmp5_1 = tmp2_0*w1;                                  EM_S[INDEX2(4,6,8)]+=tmp0 + tmp1;
7864                                  const double tmp10_1 = tmp4_0*w2;                                  EM_S[INDEX2(6,0,8)]+=tmp10;
7865                                  const double tmp11_1 = d_2*w3;                                  EM_S[INDEX2(6,2,8)]+=tmp2 + tmp3;
7866                                  const double tmp0_1 = tmp0_0*w0;                                  EM_S[INDEX2(6,4,8)]+=tmp0 + tmp1;
7867                                  EM_S[INDEX2(0,0,8)]+=tmp13_1 + tmp14_1 + tmp15_1;                                  EM_S[INDEX2(6,6,8)]+=d_0*w3 + d_3*w4 + tmp9;
7868                                  EM_S[INDEX2(2,0,8)]+=tmp6_1 + tmp7_1;                              } else { // constant data
7869                                  EM_S[INDEX2(4,0,8)]+=tmp4_1 + tmp5_1;                                  const double d_0 = d_p[0];
7870                                  EM_S[INDEX2(6,0,8)]+=tmp18_1;                                  EM_S[INDEX2(0,0,8)]+=16*d_0*w2;
7871                                  EM_S[INDEX2(0,2,8)]+=tmp6_1 + tmp7_1;                                  EM_S[INDEX2(0,2,8)]+=8*d_0*w2;
7872                                  EM_S[INDEX2(2,2,8)]+=tmp10_1 + tmp11_1 + tmp12_1;                                  EM_S[INDEX2(0,4,8)]+=8*d_0*w2;
7873                                  EM_S[INDEX2(4,2,8)]+=tmp18_1;                                  EM_S[INDEX2(0,6,8)]+=4*d_0*w2;
7874                                  EM_S[INDEX2(6,2,8)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX2(2,0,8)]+=8*d_0*w2;
7875                                  EM_S[INDEX2(0,4,8)]+=tmp4_1 + tmp5_1;                                  EM_S[INDEX2(2,2,8)]+=16*d_0*w2;
7876                                  EM_S[INDEX2(2,4,8)]+=tmp18_1;                                  EM_S[INDEX2(2,4,8)]+=4*d_0*w2;
7877                                  EM_S[INDEX2(4,4,8)]+=tmp10_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX2(2,6,8)]+=8*d_0*w2;
7878                                  EM_S[INDEX2(6,4,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(4,0,8)]+=8*d_0*w2;
7879                                  EM_S[INDEX2(0,6,8)]+=tmp18_1;                                  EM_S[INDEX2(4,2,8)]+=4*d_0*w2;
7880                                  EM_S[INDEX2(2,6,8)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX2(4,4,8)]+=16*d_0*w2;
7881                                  EM_S[INDEX2(4,6,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(4,6,8)]+=8*d_0*w2;
7882                                  EM_S[INDEX2(6,6,8)]+=tmp15_1 + tmp16_1 + tmp17_1;                                  EM_S[INDEX2(6,0,8)]+=4*d_0*w2;
7883                              } else { /* constant data */                                  EM_S[INDEX2(6,2,8)]+=8*d_0*w2;
7884                                  const double tmp0_1 = d_p[0]*w5;                                  EM_S[INDEX2(6,4,8)]+=8*d_0*w2;
7885                                  const double tmp2_1 = d_p[0]*w7;                                  EM_S[INDEX2(6,6,8)]+=16*d_0*w2;
                                 const double tmp1_1 = d_p[0]*w6;  
                                 EM_S[INDEX2(0,0,8)]+=tmp1_1;  
                                 EM_S[INDEX2(2,0,8)]+=tmp0_1;  
                                 EM_S[INDEX2(4,0,8)]+=tmp0_1;  
                                 EM_S[INDEX2(6,0,8)]+=tmp2_1;  
                                 EM_S[INDEX2(0,2,8)]+=tmp0_1;  
                                 EM_S[INDEX2(2,2,8)]+=tmp1_1;  
                                 EM_S[INDEX2(4,2,8)]+=tmp2_1;  
                                 EM_S[INDEX2(6,2,8)]+=tmp0_1;  
                                 EM_S[INDEX2(0,4,8)]+=tmp0_1;  
                                 EM_S[INDEX2(2,4,8)]+=tmp2_1;  
                                 EM_S[INDEX2(4,4,8)]+=tmp1_1;  
                                 EM_S[INDEX2(6,4,8)]+=tmp0_1;  
                                 EM_S[INDEX2(0,6,8)]+=tmp2_1;  
                                 EM_S[INDEX2(2,6,8)]+=tmp0_1;  
                                 EM_S[INDEX2(4,6,8)]+=tmp0_1;  
                                 EM_S[INDEX2(6,6,8)]+=tmp1_1;  
7886                              }                              }
7887                          }                          }
7888                          ///////////////                          ///////////////
# Line 7947  void Brick::assemblePDEBoundarySingle(Pa Line 7895  void Brick::assemblePDEBoundarySingle(Pa
7895                                  const double y_1 = y_p[1];                                  const double y_1 = y_p[1];
7896                                  const double y_2 = y_p[2];                                  const double y_2 = y_p[2];
7897                                  const double y_3 = y_p[3];                                  const double y_3 = y_p[3];
7898                                  const double tmp0_0 = y_1 + y_2;                                  const double tmp0 = 6*w2*(y_1 + y_2);
7899                                  const double tmp1_0 = y_0 + y_3;                                  const double tmp1 = 6*w2*(y_0 + y_3);
7900                                  const double tmp2_1 = w10*y_3;                                  EM_F[0]+=tmp0 + 6*w0*y_3 + 6*w1*y_0;
7901                                  const double tmp8_1 = w10*y_0;                                  EM_F[2]+=tmp1 + 6*w0*y_2 + 6*w1*y_1;
7902                                  const double tmp5_1 = w10*y_2;                                  EM_F[4]+=tmp1 + 6*w0*y_1 + 6*w1*y_2;
7903                                  const double tmp3_1 = w8*y_1;                                  EM_F[6]+=tmp0 + 6*w0*y_0 + 6*w1*y_3;
7904                                  const double tmp9_1 = w8*y_3;                              } else { // constant data
7905                                  const double tmp0_1 = w8*y_0;                                  EM_F[0]+=36*w2*y_p[0];
7906                                  const double tmp1_1 = tmp0_0*w9;                                  EM_F[2]+=36*w2*y_p[0];
7907                                  const double tmp7_1 = w8*y_2;                                  EM_F[4]+=36*w2*y_p[0];
7908                                  const double tmp4_1 = tmp1_0*w9;                                  EM_F[6]+=36*w2*y_p[0];
                                 const double tmp6_1 = w10*y_1;  
                                 EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;  
                                 EM_F[2]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                 EM_F[4]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                 EM_F[6]+=tmp1_1 + tmp8_1 + tmp9_1;  
                             } else { /* constant data */  
                                 const double tmp0_1 = w11*y_p[0];  
                                 EM_F[0]+=tmp0_1;  
                                 EM_F[2]+=tmp0_1;  
                                 EM_F[4]+=tmp0_1;  
                                 EM_F[6]+=tmp0_1;  
7909                              }                              }
7910                          }                          }
7911                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*k1;                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*k1;
# Line 7997  void Brick::assemblePDEBoundarySingle(Pa Line 7934  void Brick::assemblePDEBoundarySingle(Pa
7934                                  const double d_1 = d_p[1];                                  const double d_1 = d_p[1];
7935                                  const double d_2 = d_p[2];                                  const double d_2 = d_p[2];
7936                                  const double d_3 = d_p[3];                                  const double d_3 = d_p[3];
7937                                  const double tmp1_0 = d_1 + d_3;                                  const double tmp0 = w0*(d_0 + d_2);
7938                                  const double tmp6_0 = d_0 + d_1 + d_2 + d_3;                                  const double tmp1 = w1*(d_1 + d_3);
7939                                  const double tmp0_0 = d_0 + d_2;                                  const double tmp2 = w0*(d_2 + d_3);
7940                                  const double tmp3_0 = d_0 + d_1;                                  const double tmp3 = w1*(d_0 + d_1);
7941                                  const double tmp4_0 = d_0 + d_3;                                  const double tmp4 = w0*(d_1 + d_3);
7942                                  const double tmp2_0 = d_2 + d_3;                                  const double tmp5 = w1*(d_0 + d_2);
7943                                  const double tmp5_0 = d_1 + d_2;                                  const double tmp6 = w2*(d_0 + d_3);
7944                                  const double tmp10_1 = d_3*w4;                                  const double tmp7 = w2*(d_1 + d_2);
7945                                  const double tmp13_1 = tmp2_0*w1;                                  const double tmp8 = w0*(d_0 + d_1);
7946                                  const double tmp16_1 = d_0*w4;                                  const double tmp9 = w1*(d_2 + d_3);
7947                                  const double tmp18_1 = d_2*w4;                                  const double tmp10 = w2*(d_0 + d_1 + d_2 + d_3);
7948                                  const double tmp12_1 = tmp3_0*w0;                                  EM_S[INDEX2(1,1,8)]+=d_0*w4 + d_3*w3 + tmp7;
7949                                  const double tmp3_1 = tmp3_0*w1;                                  EM_S[INDEX2(1,3,8)]+=tmp2 + tmp3;
7950                                  const double tmp5_1 = tmp0_0*w1;                                  EM_S[INDEX2(1,5,8)]+=tmp4 + tmp5;
7951                                  const double tmp17_1 = d_1*w3;                                  EM_S[INDEX2(1,7,8)]+=tmp10;
7952                                  const double tmp9_1 = d_0*w3;                                  EM_S[INDEX2(3,1,8)]+=tmp2 + tmp3;
7953                                  const double tmp14_1 = tmp6_0*w2;                                  EM_S[INDEX2(3,3,8)]+=d_1*w4 + d_2*w3 + tmp6;
7954                                  const double tmp1_1 = tmp1_0*w1;                                  EM_S[INDEX2(3,5,8)]+=tmp10;
7955                                  const double tmp11_1 = tmp5_0*w2;                                  EM_S[INDEX2(3,7,8)]+=tmp0 + tmp1;
7956                                  const double tmp4_1 = tmp1_0*w0;                                  EM_S[INDEX2(5,1,8)]+=tmp4 + tmp5;
7957                                  const double tmp2_1 = tmp2_0*w0;                                  EM_S[INDEX2(5,3,8)]+=tmp10;
7958                                  const double tmp15_1 = d_3*w3;                                  EM_S[INDEX2(5,5,8)]+=d_1*w3 + d_2*w4 + tmp6;
7959                                  const double tmp7_1 = d_1*w4;                                  EM_S[INDEX2(5,7,8)]+=tmp8 + tmp9;
7960                                  const double tmp8_1 = tmp4_0*w2;                                  EM_S[INDEX2(7,1,8)]+=tmp10;
7961                                  const double tmp6_1 = d_2*w3;                                  EM_S[INDEX2(7,3,8)]+=tmp0 + tmp1;
7962                                  const double tmp0_1 = tmp0_0*w0;                                  EM_S[INDEX2(7,5,8)]+=tmp8 + tmp9;
7963                                  EM_S[INDEX2(1,1,8)]+=tmp11_1 + tmp15_1 + tmp16_1;                                  EM_S[INDEX2(7,7,8)]+=d_0*w3 + d_3*w4 + tmp7;
7964                                  EM_S[INDEX2(3,1,8)]+=tmp2_1 + tmp3_1;                              } else { // constant data
7965                                  EM_S[INDEX2(5,1,8)]+=tmp4_1 + tmp5_1;                                  const double d_0 = d_p[0];
7966                                  EM_S[INDEX2(7,1,8)]+=tmp14_1;                                  EM_S[INDEX2(1,1,8)]+=16*d_0*w2;
7967                                  EM_S[INDEX2(1,3,8)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX2(1,3,8)]+=8*d_0*w2;
7968                                  EM_S[INDEX2(3,3,8)]+=tmp6_1 + tmp7_1 + tmp8_1;                                  EM_S[INDEX2(1,5,8)]+=8*d_0*w2;
7969                                  EM_S[INDEX2(5,3,8)]+=tmp14_1;                                  EM_S[INDEX2(1,7,8)]+=4*d_0*w2;
7970                                  EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(3,1,8)]+=8*d_0*w2;
7971                                  EM_S[INDEX2(1,5,8)]+=tmp4_1 + tmp5_1;                                  EM_S[INDEX2(3,3,8)]+=16*d_0*w2;
7972                                  EM_S[INDEX2(3,5,8)]+=tmp14_1;                                  EM_S[INDEX2(3,5,8)]+=4*d_0*w2;
7973                                  EM_S[INDEX2(5,5,8)]+=tmp17_1 + tmp18_1 + tmp8_1;                                  EM_S[INDEX2(3,7,8)]+=8*d_0*w2;
7974                                  EM_S[INDEX2(7,5,8)]+=tmp12_1 + tmp13_1;                                  EM_S[INDEX2(5,1,8)]+=8*d_0*w2;
7975                                  EM_S[INDEX2(1,7,8)]+=tmp14_1;                                  EM_S[INDEX2(5,3,8)]+=4*d_0*w2;
7976                                  EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(5,5,8)]+=16*d_0*w2;
7977                                  EM_S[INDEX2(5,7,8)]+=tmp12_1 + tmp13_1;                                  EM_S[INDEX2(5,7,8)]+=8*d_0*w2;
7978                                  EM_S[INDEX2(7,7,8)]+=tmp10_1 + tmp11_1 + tmp9_1;                                  EM_S[INDEX2(7,1,8)]+=4*d_0*w2;
7979                              } else { /* constant data */                                  EM_S[INDEX2(7,3,8)]+=8*d_0*w2;
7980                                  const double tmp0_1 = d_p[0]*w5;                                  EM_S[INDEX2(7,5,8)]+=8*d_0*w2;
7981                                  const double tmp2_1 = d_p[0]*w7;                                  EM_S[INDEX2(7,7,8)]+=16*d_0*w2;
                                 const double tmp1_1 = d_p[0]*w6;  
                                 EM_S[INDEX2(1,1,8)]+=tmp1_1;  
                                 EM_S[INDEX2(3,1,8)]+=tmp0_1;  
                                 EM_S[INDEX2(5,1,8)]+=tmp0_1;  
                                 EM_S[INDEX2(7,1,8)]+=tmp2_1;  
                                 EM_S[INDEX2(1,3,8)]+=tmp0_1;  
                                 EM_S[INDEX2(3,3,8)]+=tmp1_1;  
                                 EM_S[INDEX2(5,3,8)]+=tmp2_1;  
                                 EM_S[INDEX2(7,3,8)]+=tmp0_1;  
                                 EM_S[INDEX2(1,5,8)]+=tmp0_1;  
                                 EM_S[INDEX2(3,5,8)]+=tmp2_1;  
                                 EM_S[INDEX2(5,5,8)]+=tmp1_1;  
                                 EM_S[INDEX2(7,5,8)]+=tmp0_1;  
                                 EM_S[INDEX2(1,7,8)]+=tmp2_1;  
                                 EM_S[INDEX2(3,7,8)]+=tmp0_1;  
                                 EM_S[INDEX2(5,7,8)]+=tmp0_1;  
                                 EM_S[INDEX2(7,7,8)]+=tmp1_1;  
7982                              }                              }
7983                          }                          }
7984                          ///////////////                          ///////////////
# Line 8071  void Brick::assemblePDEBoundarySingle(Pa Line 7991  void Brick::assemblePDEBoundarySingle(Pa
7991                                  const double y_1 = y_p[1];                                  const double y_1 = y_p[1];
7992                                  const double y_2 = y_p[2];                                  const double y_2 = y_p[2];
7993                                  const double y_3 = y_p[3];                                  const double y_3 = y_p[3];
7994                                  const double tmp0_0 = y_1 + y_2;                                  const double tmp0 = 6*w2*(y_1 + y_2);
7995                                  const double tmp1_0 = y_0 + y_3;                                  const double tmp1 = 6*w2*(y_0 + y_3);
7996                                  const double tmp2_1 = w10*y_3;                                  EM_F[1]+=tmp0 + 6*w0*y_3 + 6*w1*y_0;
7997                                  const double tmp8_1 = w10*y_0;                                  EM_F[3]+=tmp1 + 6*w0*y_2 + 6*w1*y_1;
7998                                  const double tmp5_1 = w10*y_2;                                  EM_F[5]+=tmp1 + 6*w0*y_1 + 6*w1*y_2;
7999                                  const double tmp3_1 = w8*y_1;                                  EM_F[7]+=tmp0 + 6*w0*y_0 + 6*w1*y_3;
8000                                  const double tmp9_1 = w8*y_3;                              } else { // constant data
8001                                  const double tmp0_1 = w8*y_0;                                  EM_F[1]+=36*w2*y_p[0];
8002                                  const double tmp1_1 = tmp0_0*w9;                                  EM_F[3]+=36*w2*y_p[0];
8003                                  const double tmp7_1 = w8*y_2;                                  EM_F[5]+=36*w2*y_p[0];
8004                                  const double tmp4_1 = tmp1_0*w9;                                  EM_F[7]+=36*w2*y_p[0];
                                 const double tmp6_1 = w10*y_1;  
                                 EM_F[1]+=tmp0_1 + tmp1_1 + tmp2_1;  
                                 EM_F[3]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                 EM_F[5]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                 EM_F[7]+=tmp1_1 + tmp8_1 + tmp9_1;  
                             } else { /* constant data */  
                                 const double tmp0_1 = w11*y_p[0];  
                                 EM_F[1]+=tmp0_1;  
                                 EM_F[3]+=tmp0_1;  
                                 EM_F[5]+=tmp0_1;  
                                 EM_F[7]+=tmp0_1;  
8005                              }                              }
8006                          }                          }
8007                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*(k1+1)-2;                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*(k1+1)-2;
# Line 8121  void Brick::assemblePDEBoundarySingle(Pa Line 8030  void Brick::assemblePDEBoundarySingle(Pa
8030                                  const double d_1 = d_p[1];                                  const double d_1 = d_p[1];
8031                                  const double d_2 = d_p[2];                                  const double d_2 = d_p[2];
8032                                  const double d_3 = d_p[3];                                  const double d_3 = d_p[3];
8033                                  const double tmp2_0 = d_1 + d_3;                                  const double tmp0 = w5*(d_0 + d_1);
8034                                  const double tmp5_0 = d_0 + d_1 + d_2 + d_3;                                  const double tmp1 = w6*(d_2 + d_3);
8035                                  const double tmp3_0 = d_0 + d_2;                                  const double tmp2 = w5*(d_0 + d_2);
8036                                  const double tmp1_0 = d_0 + d_1;                                  const double tmp3 = w6*(d_1 + d_3);
8037                                  const double tmp4_0 = d_0 + d_3;                                  const double tmp4 = w5*(d_1 + d_3);
8038                                  const double tmp0_0 = d_2 + d_3;                                  const double tmp5 = w6*(d_0 + d_2);
8039                                  const double tmp6_0 = d_1 + d_2;                                  const double tmp6 = w7*(d_0 + d_3);
8040                                  const double tmp2_1 = tmp2_0*w13;                                  const double tmp7 = w7*(d_0 + d_1 + d_2 + d_3);
8041                                  const double tmp14_1 = d_3*w15;                                  const double tmp8 = w7*(d_1 + d_2);
8042                                  const double tmp0_1 = tmp0_0*w13;                                  const double tmp9 = w5*(d_2 + d_3);
8043                                  const double tmp3_1 = tmp3_0*w12;                                  const double tmp10 = w6*(d_0 + d_1);
8044                                  const double tmp17_1 = tmp1_0*w13;                                  EM_S[INDEX2(0,0,8)]+=d_0*w9 + d_3*w8 + tmp8;
8045                                  const double tmp18_1 = tmp0_0*w12;                                  EM_S[INDEX2(0,1,8)]+=tmp10 + tmp9;
8046                                  const double tmp8_1 = d_1*w15;                                  EM_S[INDEX2(0,4,8)]+=tmp4 + tmp5;
8047                                  const double tmp16_1 = d_0*w15;                                  EM_S[INDEX2(0,5,8)]+=tmp7;
8048                                  const double tmp11_1 = d_2*w15;                                  EM_S[INDEX2(1,0,8)]+=tmp10 + tmp9;
8049                                  const double tmp5_1 = tmp2_0*w12;                                  EM_S[INDEX2(1,1,8)]+=d_1*w9 + d_2*w8 + tmp6;
8050                                  const double tmp15_1 = d_3*w16;                                  EM_S[INDEX2(1,4,8)]+=tmp7;
8051                                  const double tmp13_1 = tmp6_0*w14;                                  EM_S[INDEX2(1,5,8)]+=tmp2 + tmp3;
8052                                  const double tmp1_1 = tmp1_0*w12;                                  EM_S[INDEX2(4,0,8)]+=tmp4 + tmp5;
8053                                  const double tmp7_1 = tmp4_0*w14;                                  EM_S[INDEX2(4,1,8)]+=tmp7;
8054                                  const double tmp12_1 = d_0*w16;                                  EM_S[INDEX2(4,4,8)]+=d_1*w8 + d_2*w9 + tmp6;
8055                                  const double tmp10_1 = d_1*w16;                                  EM_S[INDEX2(4,5,8)]+=tmp0 + tmp1;
8056                                  const double tmp6_1 = d_2*w16;                                  EM_S[INDEX2(5,0,8)]+=tmp7;
8057                                  const double tmp9_1 = tmp5_0*w14;                                  EM_S[INDEX2(5,1,8)]+=tmp2 + tmp3;
8058                                  const double tmp4_1 = tmp3_0*w13;                                  EM_S[INDEX2(5,4,8)]+=tmp0 + tmp1;
8059                                  EM_S[INDEX2(0,0,8)]+=tmp12_1 + tmp13_1 + tmp14_1;                                  EM_S[INDEX2(5,5,8)]+=d_0*w8 + d_3*w9 + tmp8;
8060                                  EM_S[INDEX2(1,0,8)]+=tmp17_1 + tmp18_1;                              } else { // constant data
8061                                  EM_S[INDEX2(4,0,8)]+=tmp4_1 + tmp5_1;                                  const double d_0 = d_p[0];
8062                                  EM_S[INDEX2(5,0,8)]+=tmp9_1;                                  EM_S[INDEX2(0,0,8)]+=16*d_0*w7;
8063                                  EM_S[INDEX2(0,1,8)]+=tmp17_1 + tmp18_1;                                  EM_S[INDEX2(0,1,8)]+=8*d_0*w7;
8064                                  EM_S[INDEX2(1,1,8)]+=tmp10_1 + tmp11_1 + tmp7_1;                                  EM_S[INDEX2(0,4,8)]+=8*d_0*w7;
8065                                  EM_S[INDEX2(4,1,8)]+=tmp9_1;                                  EM_S[INDEX2(0,5,8)]+=4*d_0*w7;
8066                                  EM_S[INDEX2(5,1,8)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX2(1,0,8)]+=8*d_0*w7;
8067                                  EM_S[INDEX2(0,4,8)]+=tmp4_1 + tmp5_1;                                  EM_S[INDEX2(1,1,8)]+=16*d_0*w7;
8068                                  EM_S[INDEX2(1,4,8)]+=tmp9_1;                                  EM_S[INDEX2(1,4,8)]+=4*d_0*w7;
8069                                  EM_S[INDEX2(4,4,8)]+=tmp6_1 + tmp7_1 + tmp8_1;                                  EM_S[INDEX2(1,5,8)]+=8*d_0*w7;
8070                                  EM_S[INDEX2(5,4,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(4,0,8)]+=8*d_0*w7;
8071                                  EM_S[INDEX2(0,5,8)]+=tmp9_1;                                  EM_S[INDEX2(4,1,8)]+=4*d_0*w7;
8072                                  EM_S[INDEX2(1,5,8)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX2(4,4,8)]+=16*d_0*w7;
8073                                  EM_S[INDEX2(4,5,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(4,5,8)]+=8*d_0*w7;
8074                                  EM_S[INDEX2(5,5,8)]+=tmp13_1 + tmp15_1 + tmp16_1;                                  EM_S[INDEX2(5,0,8)]+=4*d_0*w7;
8075                              } else { /* constant data */                                  EM_S[INDEX2(5,1,8)]+=8*d_0*w7;
8076                                  const double tmp0_1 = d_p[0]*w17;                                  EM_S[INDEX2(5,4,8)]+=8*d_0*w7;
8077                                  const double tmp2_1 = d_p[0]*w19;                                  EM_S[INDEX2(5,5,8)]+=16*d_0*w7;
                                 const double tmp1_1 = d_p[0]*w18;  
                                 EM_S[INDEX2(0,0,8)]+=tmp1_1;  
                                 EM_S[INDEX2(1,0,8)]+=tmp0_1;  
                                 EM_S[INDEX2(4,0,8)]+=tmp0_1;  
                                 EM_S[INDEX2(5,0,8)]+=tmp2_1;  
                                 EM_S[INDEX2(0,1,8)]+=tmp0_1;  
                                 EM_S[INDEX2(1,1,8)]+=tmp1_1;  
                                 EM_S[INDEX2(4,1,8)]+=tmp2_1;  
                                 EM_S[INDEX2(5,1,8)]+=tmp0_1;  
                                 EM_S[INDEX2(0,4,8)]+=tmp0_1;  
                                 EM_S[INDEX2(1,4,8)]+=tmp2_1;  
                                 EM_S[INDEX2(4,4,8)]+=tmp1_1;  
                                 EM_S[INDEX2(5,4,8)]+=tmp0_1;  
                                 EM_S[INDEX2(0,5,8)]+=tmp2_1;  
                                 EM_S[INDEX2(1,5,8)]+=tmp0_1;  
                                 EM_S[INDEX2(4,5,8)]+=tmp0_1;  
                                 EM_S[INDEX2(5,5,8)]+=tmp1_1;  
8078                              }                              }
8079                          }                          }
8080                          ///////////////                          ///////////////
# Line 8195  void Brick::assemblePDEBoundarySingle(Pa Line 8087  void Brick::assemblePDEBoundarySingle(Pa
8087                                  const double y_1 = y_p[1];                                  const double y_1 = y_p[1];
8088                                  const double y_2 = y_p[2];                                  const double y_2 = y_p[2];
8089                                  const double y_3 = y_p[3];                                  const double y_3 = y_p[3];
8090                                  const double tmp0_0 = y_1 + y_2;                                  const double tmp0 = 6*w7*(y_1 + y_2);
8091                                  const double tmp1_0 = y_0 + y_3;                                  const double tmp1 = 6*w7*(y_0 + y_3);
8092                                  const double tmp0_1 = w22*y_3;                                  EM_F[0]+=tmp0 + 6*w5*y_3 + 6*w6*y_0;
8093                                  const double tmp6_1 = w22*y_1;                                  EM_F[1]+=tmp1 + 6*w5*y_2 + 6*w6*y_1;
8094                                  const double tmp3_1 = w22*y_2;                                  EM_F[4]+=tmp1 + 6*w5*y_1 + 6*w6*y_2;
8095                                  const double tmp5_1 = w20*y_1;                                  EM_F[5]+=tmp0 + 6*w5*y_0 + 6*w6*y_3;
8096                                  const double tmp9_1 = w20*y_3;                              } else { // constant data
8097                                  const double tmp4_1 = tmp1_0*w21;                                  EM_F[0]+=36*w7*y_p[0];
8098                                  const double tmp8_1 = w22*y_0;                                  EM_F[1]+=36*w7*y_p[0];
8099                                  const double tmp2_1 = w20*y_0;                                  EM_F[4]+=36*w7*y_p[0];
8100                                  const double tmp7_1 = w20*y_2;                                  EM_F[5]+=36*w7*y_p[0];
                                 const double tmp1_1 = tmp0_0*w21;  
                                 EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;  
                                 EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                 EM_F[4]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                 EM_F[5]+=tmp1_1 + tmp8_1 + tmp9_1;  
                             } else { /* constant data */  
                                 const double tmp0_1 = w23*y_p[0];  
                                 EM_F[0]+=tmp0_1;  
                                 EM_F[1]+=tmp0_1;  
                                 EM_F[4]+=tmp0_1;  
                                 EM_F[5]+=tmp0_1;  
8101                              }                              }
8102                          }                          }
8103                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+k0;                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+k0;
# Line 8245  void Brick::assemblePDEBoundarySingle(Pa Line 8126  void Brick::assemblePDEBoundarySingle(Pa
8126                                  const double d_1 = d_p[1];                                  const double d_1 = d_p[1];
8127                                  const double d_2 = d_p[2];                                  const double d_2 = d_p[2];
8128                                  const double d_3 = d_p[3];                                  const double d_3 = d_p[3];
8129                                  const double tmp0_0 = d_1 + d_3;                                  const double tmp0 = w5*(d_0 + d_2);
8130                                  const double tmp2_0 = d_0 + d_1 + d_2 + d_3;                                  const double tmp1 = w6*(d_1 + d_3);
8131                                  const double tmp1_0 = d_0 + d_2;                                  const double tmp2 = w5*(d_1 + d_3);
8132                                  const double tmp4_0 = d_0 + d_1;                                  const double tmp3 = w6*(d_0 + d_2);
8133                                  const double tmp5_0 = d_0 + d_3;                                  const double tmp4 = w7*(d_0 + d_1 + d_2 + d_3);
8134                                  const double tmp3_0 = d_2 + d_3;                                  const double tmp5 = w5*(d_0 + d_1);
8135                                  const double tmp6_0 = d_1 + d_2;                                  const double tmp6 = w6*(d_2 + d_3);
8136                                  const double tmp15_1 = tmp4_0*w13;                                  const double tmp7 = w7*(d_0 + d_3);
8137                                  const double tmp10_1 = d_0*w16;                                  const double tmp8 = w7*(d_1 + d_2);
8138                                  const double tmp6_1 = tmp4_0*w12;                                  const double tmp9 = w5*(d_2 + d_3);
8139                                  const double tmp16_1 = tmp3_0*w12;                                  const double tmp10 = w6*(d_0 + d_1);
8140                                  const double tmp0_1 = tmp0_0*w13;                                  EM_S[INDEX2(2,2,8)]+=d_0*w9 + d_3*w8 + tmp8;
8141                                  const double tmp2_1 = tmp1_0*w13;                                  EM_S[INDEX2(2,3,8)]+=tmp10 + tmp9;
8142                                  const double tmp18_1 = d_1*w15;                                  EM_S[INDEX2(2,6,8)]+=tmp2 + tmp3;
8143                                  const double tmp14_1 = d_0*w15;                                  EM_S[INDEX2(2,7,8)]+=tmp4;
8144                                  const double tmp9_1 = d_2*w15;                                  EM_S[INDEX2(3,2,8)]+=tmp10 + tmp9;
8145                                  const double tmp4_1 = tmp2_0*w14;                                  EM_S[INDEX2(3,3,8)]+=d_1*w9 + d_2*w8 + tmp7;
8146                                  const double tmp13_1 = d_3*w16;                                  EM_S[INDEX2(3,6,8)]+=tmp4;
8147                                  const double tmp11_1 = tmp6_0*w14;                                  EM_S[INDEX2(3,7,8)]+=tmp0 + tmp1;
8148                                  const double tmp1_1 = tmp1_0*w12;                                  EM_S[INDEX2(6,2,8)]+=tmp2 + tmp3;
8149                                  const double tmp12_1 = d_3*w15;                                  EM_S[INDEX2(6,3,8)]+=tmp4;
8150                                  const double tmp3_1 = tmp0_0*w12;                                  EM_S[INDEX2(6,6,8)]+=d_1*w8 + d_2*w9 + tmp7;
8151                                  const double tmp7_1 = d_1*w16;                                  EM_S[INDEX2(6,7,8)]+=tmp5 + tmp6;
8152                                  const double tmp17_1 = d_2*w16;                                  EM_S[INDEX2(7,2,8)]+=tmp4;
8153                                  const double tmp8_1 = tmp5_0*w14;                                  EM_S[INDEX2(7,3,8)]+=tmp0 + tmp1;
8154                                  const double tmp5_1 = tmp3_0*w13;                                  EM_S[INDEX2(7,6,8)]+=tmp5 + tmp6;
8155                                  EM_S[INDEX2(2,2,8)]+=tmp10_1 + tmp11_1 + tmp12_1;                                  EM_S[INDEX2(7,7,8)]+=d_0*w8 + d_3*w9 + tmp8;
8156                                  EM_S[INDEX2(3,2,8)]+=tmp15_1 + tmp16_1;                              } else { // constant data
8157                                  EM_S[INDEX2(6,2,8)]+=tmp2_1 + tmp3_1;                                  const double d_0 = d_p[0];
8158                                  EM_S[INDEX2(7,2,8)]+=tmp4_1;                                  EM_S[INDEX2(2,2,8)]+=16*d_0*w7;
8159                                  EM_S[INDEX2(2,3,8)]+=tmp15_1 + tmp16_1;                                  EM_S[INDEX2(2,3,8)]+=8*d_0*w7;
8160                                  EM_S[INDEX2(3,3,8)]+=tmp7_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX2(2,6,8)]+=8*d_0*w7;
8161                                  EM_S[INDEX2(6,3,8)]+=tmp4_1;                                  EM_S[INDEX2(2,7,8)]+=4*d_0*w7;
8162                                  EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(3,2,8)]+=8*d_0*w7;
8163                                  EM_S[INDEX2(2,6,8)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX2(3,3,8)]+=16*d_0*w7;
8164                                  EM_S[INDEX2(3,6,8)]+=tmp4_1;                                  EM_S[INDEX2(3,6,8)]+=4*d_0*w7;
8165                                  EM_S[INDEX2(6,6,8)]+=tmp17_1 + tmp18_1 + tmp8_1;                                  EM_S[INDEX2(3,7,8)]+=8*d_0*w7;
8166                                  EM_S[INDEX2(7,6,8)]+=tmp5_1 + tmp6_1;                                  EM_S[INDEX2(6,2,8)]+=8*d_0*w7;
8167                                  EM_S[INDEX2(2,7,8)]+=tmp4_1;                                  EM_S[INDEX2(6,3,8)]+=4*d_0*w7;
8168                                  EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(6,6,8)]+=16*d_0*w7;
8169                                  EM_S[INDEX2(6,7,8)]+=tmp5_1 + tmp6_1;                                  EM_S[INDEX2(6,7,8)]+=8*d_0*w7;
8170                                  EM_S[INDEX2(7,7,8)]+=tmp11_1 + tmp13_1 + tmp14_1;                                  EM_S[INDEX2(7,2,8)]+=4*d_0*w7;
8171                              } else { /* constant data */                                  EM_S[INDEX2(7,3,8)]+=8*d_0*w7;
8172                                  const double tmp0_1 = d_p[0]*w17;                                  EM_S[INDEX2(7,6,8)]+=8*d_0*w7;
8173                                  const double tmp1_1 = d_p[0]*w19;                                  EM_S[INDEX2(7,7,8)]+=16*d_0*w7;
                                 const double tmp2_1 = d_p[0]*w18;  
                                 EM_S[INDEX2(2,2,8)]+=tmp2_1;  
                                 EM_S[INDEX2(3,2,8)]+=tmp0_1;  
                                 EM_S[INDEX2(6,2,8)]+=tmp0_1;  
                                 EM_S[INDEX2(7,2,8)]+=tmp1_1;  
                                 EM_S[INDEX2(2,3,8)]+=tmp0_1;  
                                 EM_S[INDEX2(3,3,8)]+=tmp2_1;  
                                 EM_S[INDEX2(6,3,8)]+=tmp1_1;  
                                 EM_S[INDEX2(7,3,8)]+=tmp0_1;  
                                 EM_S[INDEX2(2,6,8)]+=tmp0_1;  
                                 EM_S[INDEX2(3,6,8)]+=tmp1_1;  
                                 EM_S[INDEX2(6,6,8)]+=tmp2_1;  
                                 EM_S[INDEX2(7,6,8)]+=tmp0_1;  
                                 EM_S[INDEX2(2,7,8)]+=tmp1_1;  
                                 EM_S[INDEX2(3,7,8)]+=tmp0_1;  
                                 EM_S[INDEX2(6,7,8)]+=tmp0_1;  
                                 EM_S[INDEX2(7,7,8)]+=tmp2_1;  
8174                              }                              }
8175                          }                          }
8176                          ///////////////                          ///////////////
# Line 8319  void Brick::assemblePDEBoundarySingle(Pa Line 8183  void Brick::assemblePDEBoundarySingle(Pa
8183                                  const double y_1 = y_p[1];                                  const double y_1 = y_p[1];
8184                                  const double y_2 = y_p[2];                                  const double y_2 = y_p[2];
8185                                  const double y_3 = y_p[3];                                  const double y_3 = y_p[3];
8186                                  const double tmp0_0 = y_1 + y_2;                                  const double tmp0 = 6*w7*(y_1 + y_2);
8187                                  const double tmp1_0 = y_0 + y_3;                                  const double tmp1 = 6*w7*(y_0 + y_3);
8188                                  const double tmp0_1 = w22*y_3;                                  EM_F[2]+=tmp0 + 6*w5*y_3 + 6*w6*y_0;
8189                                  const double tmp6_1 = w22*y_1;                                  EM_F[3]+=tmp1 + 6*w5*y_2 + 6*w6*y_1;
8190                                  const double tmp3_1 = w22*y_2;                                  EM_F[6]+=tmp1 + 6*w5*y_1 + 6*w6*y_2;
8191                                  const double tmp5_1 = w20*y_1;                                  EM_F[7]+=tmp0 + 6*w5*y_0 + 6*w6*y_3;
8192                                  const double tmp9_1 = w20*y_3;                              } else { // constant data
8193                                  const double tmp4_1 = tmp1_0*w21;                                  EM_F[2]+=36*w7*y_p[0];
8194                                  const double tmp8_1 = w22*y_0;                                  EM_F[3]+=36*w7*y_p[0];
8195                                  const double tmp2_1 = w20*y_0;                                  EM_F[6]+=36*w7*y_p[0];
8196                                  const double tmp7_1 = w20*y_2;                                  EM_F[7]+=36*w7*y_p[0];
                                 const double tmp1_1 = tmp0_0*w21;  
                                 EM_F[2]+=tmp0_1 + tmp1_1 + tmp2_1;  
                                 EM_F[3]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                 EM_F[6]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                 EM_F[7]+=tmp1_1 + tmp8_1 + tmp9_1;  
                             } else { /* constant data */  
                                 const double tmp0_1 = w23*y_p[0];  
                                 EM_F[2]+=tmp0_1;  
                                 EM_F[3]+=tmp0_1;  
                                 EM_F[6]+=tmp0_1;  
                                 EM_F[7]+=tmp0_1;  
8197                              }                              }
8198                          }                          }
8199                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*(m_NN[1]-2)+k0;                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*(m_NN[1]-2)+k0;
# Line 8369  void Brick::assemblePDEBoundarySingle(Pa Line 8222  void Brick::assemblePDEBoundarySingle(Pa
8222                                  const double d_1 = d_p[1];                                  const double d_1 = d_p[1];
8223                                  const double d_2 = d_p[2];                                  const double d_2 = d_p[2];
8224                                  const double d_3 = d_p[3];                                  const double d_3 = d_p[3];
8225                                  const double tmp0_0 = d_1 + d_3;                                  const double tmp0 = w10*(d_0 + d_2);
8226                                  const double tmp2_0 = d_0 + d_1 + d_2 + d_3;                                  const double tmp1 = w11*(d_1 + d_3);
8227                                  const double tmp1_0 = d_0 + d_2;                                  const double tmp2 = w12*(d_0 + d_1 + d_2 + d_3);
8228                                  const double tmp6_0 = d_0 + d_1;                                  const double tmp3 = w12*(d_1 + d_2);
8229                                  const double tmp4_0 = d_0 + d_3;                                  const double tmp4 = w10*(d_1 + d_3);
8230                                  const double tmp5_0 = d_2 + d_3;                                  const double tmp5 = w11*(d_0 + d_2);
8231                                  const double tmp3_0 = d_1 + d_2;                                  const double tmp6 = w12*(d_0 + d_3);
8232                                  const double tmp18_1 = tmp5_0*w24;                                  const double tmp7 = w10*(d_0 + d_1);
8233                                  const double tmp6_1 = tmp1_0*w25;                                  const double tmp8 = w11*(d_2 + d_3);
8234                                  const double tmp4_1 = d_0*w27;                                  const double tmp9 = w10*(d_2 + d_3);
8235                                  const double tmp12_1 = d_2*w27;                                  const double tmp10 = w11*(d_0 + d_1);
8236                                  const double tmp0_1 = tmp0_0*w25;                                  EM_S[INDEX2(0,0,8)]+=d_0*w14 + d_3*w13 + tmp3;
8237                                  const double tmp5_1 = tmp3_0*w26;                                  EM_S[INDEX2(0,1,8)]+=tmp10 + tmp9;
8238                                  const double tmp2_1 = tmp2_0*w26;                                  EM_S[INDEX2(0,2,8)]+=tmp4 + tmp5;
8239                                  const double tmp17_1 = tmp6_0*w25;                                  EM_S[INDEX2(0,3,8)]+=tmp2;
8240                                  const double tmp14_1 = tmp6_0*w24;                                  EM_S[INDEX2(1,0,8)]+=tmp10 + tmp9;
8241                                  const double tmp11_1 = d_1*w28;                                  EM_S[INDEX2(1,1,8)]+=d_1*w14 + d_2*w13 + tmp6;
8242                                  const double tmp9_1 = d_1*w27;                                  EM_S[INDEX2(1,2,8)]+=tmp2;
8243                                  const double tmp16_1 = d_3*w27;                                  EM_S[INDEX2(1,3,8)]+=tmp0 + tmp1;
8244                                  const double tmp8_1 = d_2*w28;                                  EM_S[INDEX2(2,0,8)]+=tmp4 + tmp5;
8245                                  const double tmp7_1 = tmp0_0*w24;                                  EM_S[INDEX2(2,1,8)]+=tmp2;
8246                                  const double tmp15_1 = d_0*w28;                                  EM_S[INDEX2(2,2,8)]+=d_1*w13 + d_2*w14 + tmp6;
8247                                  const double tmp13_1 = tmp5_0*w25;                                  EM_S[INDEX2(2,3,8)]+=tmp7 + tmp8;
8248                                  const double tmp3_1 = d_3*w28;                                  EM_S[INDEX2(3,0,8)]+=tmp2;
8249                                  const double tmp10_1 = tmp4_0*w26;                                  EM_S[INDEX2(3,1,8)]+=tmp0 + tmp1;
8250                                  const double tmp1_1 = tmp1_0*w24;                                  EM_S[INDEX2(3,2,8)]+=tmp7 + tmp8;
8251                                  EM_S[INDEX2(0,0,8)]+=tmp15_1 + tmp16_1 + tmp5_1;                                  EM_S[INDEX2(3,3,8)]+=d_0*w13 + d_3*w14 + tmp3;
8252                                  EM_S[INDEX2(1,0,8)]+=tmp17_1 + tmp18_1;                              } else { // constant data
8253                                  EM_S[INDEX2(2,0,8)]+=tmp6_1 + tmp7_1;                                  const double d_0 = d_p[0];
8254                                  EM_S[INDEX2(3,0,8)]+=tmp2_1;                                  EM_S[INDEX2(0,0,8)]+=16*d_0*w12;
8255                                  EM_S[INDEX2(0,1,8)]+=tmp17_1 + tmp18_1;                                  EM_S[INDEX2(0,1,8)]+=8*d_0*w12;
8256                                  EM_S[INDEX2(1,1,8)]+=tmp10_1 + tmp11_1 + tmp12_1;                                  EM_S[INDEX2(0,2,8)]+=8*d_0*w12;
8257                                  EM_S[INDEX2(2,1,8)]+=tmp2_1;                                  EM_S[INDEX2(0,3,8)]+=4*d_0*w12;
8258                                  EM_S[INDEX2(3,1,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(1,0,8)]+=8*d_0*w12;
8259                                  EM_S[INDEX2(0,2,8)]+=tmp6_1 + tmp7_1;                                  EM_S[INDEX2(1,1,8)]+=16*d_0*w12;
8260                                  EM_S[INDEX2(1,2,8)]+=tmp2_1;                                  EM_S[INDEX2(1,2,8)]+=4*d_0*w12;
8261                                  EM_S[INDEX2(2,2,8)]+=tmp10_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX2(1,3,8)]+=8*d_0*w12;
8262                                  EM_S[INDEX2(3,2,8)]+=tmp13_1 + tmp14_1;                                  EM_S[INDEX2(2,0,8)]+=8*d_0*w12;
8263                                  EM_S[INDEX2(0,3,8)]+=tmp2_1;                                  EM_S[INDEX2(2,1,8)]+=4*d_0*w12;
8264                                  EM_S[INDEX2(1,3,8)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX2(2,2,8)]+=16*d_0*w12;
8265                                  EM_S[INDEX2(2,3,8)]+=tmp13_1 + tmp14_1;                                  EM_S[INDEX2(2,3,8)]+=8*d_0*w12;
8266                                  EM_S[INDEX2(3,3,8)]+=tmp3_1 + tmp4_1 + tmp5_1;                                  EM_S[INDEX2(3,0,8)]+=4*d_0*w12;
8267                              } else { /* constant data */                                  EM_S[INDEX2(3,1,8)]+=8*d_0*w12;
8268                                  const double tmp2_1 = d_p[0]*w31;                                  EM_S[INDEX2(3,2,8)]+=8*d_0*w12;
8269                                  const double tmp1_1 = d_p[0]*w30;                                  EM_S[INDEX2(3,3,8)]+=16*d_0*w12;
                                 const double tmp0_1 = d_p[0]*w29;  
                                 EM_S[INDEX2(0,0,8)]+=tmp2_1;  
                                 EM_S[INDEX2(1,0,8)]+=tmp0_1;  
                                 EM_S[INDEX2(2,0,8)]+=tmp0_1;  
                                 EM_S[INDEX2(3,0,8)]+=tmp1_1;  
                                 EM_S[INDEX2(0,1,8)]+=tmp0_1;  
                                 EM_S[INDEX2(1,1,8)]+=tmp2_1;  
                                 EM_S[INDEX2(2,1,8)]+=tmp1_1;  
                                 EM_S[INDEX2(3,1,8)]+=tmp0_1;  
                                 EM_S[INDEX2(0,2,8)]+=tmp0_1;  
                                 EM_S[INDEX2(1,2,8)]+=tmp1_1;  
                                 EM_S[INDEX2(2,2,8)]+=tmp2_1;  
                                 EM_S[INDEX2(3,2,8)]+=tmp0_1;  
                                 EM_S[INDEX2(0,3,8)]+=tmp1_1;  
                                 EM_S[INDEX2(1,3,8)]+=tmp0_1;  
                                 EM_S[INDEX2(2,3,8)]+=tmp0_1;  
                                 EM_S[INDEX2(3,3,8)]+=tmp2_1;  
8270                              }                              }
8271                          }                          }
8272                          ///////////////                          ///////////////
# Line 8443  void Brick::assemblePDEBoundarySingle(Pa Line 8279  void Brick::assemblePDEBoundarySingle(Pa
8279                                  const double y_1 = y_p[1];                                  const double y_1 = y_p[1];
8280                                  const double y_2 = y_p[2];                                  const double y_2 = y_p[2];
8281                                  const double y_3 = y_p[3];                                  const double y_3 = y_p[3];
8282                                  const double tmp0_0 = y_1 + y_2;                                  const double tmp0 = 6*w12*(y_1 + y_2);
8283                                  const double tmp1_0 = y_0 + y_3;                                  const double tmp1 = 6*w12*(y_0 + y_3);
8284                                  const double tmp7_1 = w34*y_1;                                  EM_F[0]+=tmp0 + 6*w10*y_3 + 6*w11*y_0;
8285                                  const double tmp3_1 = w32*y_1;                                  EM_F[1]+=tmp1 + 6*w10*y_2 + 6*w11*y_1;
8286                                  const double tmp8_1 = w32*y_3;                                  EM_F[2]+=tmp1 + 6*w10*y_1 + 6*w11*y_2;
8287                                  const double tmp0_1 = w32*y_0;                                  EM_F[3]+=tmp0 + 6*w10*y_0 + 6*w11*y_3;
8288                                  const double tmp2_1 = w34*y_3;                              } else { // constant data
8289                                  const double tmp9_1 = w34*y_0;                                  EM_F[0]+=36*w12*y_p[0];
8290                                  const double tmp6_1 = w32*y_2;                                  EM_F[1]+=36*w12*y_p[0];
8291                                  const double tmp5_1 = w34*y_2;                                  EM_F[2]+=36*w12*y_p[0];
8292                                  const double tmp1_1 = tmp0_0*w33;                                  EM_F[3]+=36*w12*y_p[0];
                                 const double tmp4_1 = tmp1_0*w33;  
                                 EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;  
                                 EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                 EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                 EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;  
                             } else { /* constant data */  
                                 const double tmp0_1 = w35*y_p[0];  
                                 EM_F[0]+=tmp0_1;  
                                 EM_F[1]+=tmp0_1;  
                                 EM_F[2]+=tmp0_1;  
                                 EM_F[3]+=tmp0_1;  
8293                              }                              }
8294                          }                          }
8295                          const index_t firstNode=m_NN[0]*k1+k0;                          const index_t firstNode=m_NN[0]*k1+k0;
# Line 8493  void Brick::assemblePDEBoundarySingle(Pa Line 8318  void Brick::assemblePDEBoundarySingle(Pa
8318                                  const double d_1 = d_p[1];                                  const double d_1 = d_p[1];
8319                                  const double d_2 = d_p[2];                                  const double d_2 = d_p[2];
8320                                  const double d_3 = d_p[3];                                  const double d_3 = d_p[3];
8321                                  const double tmp2_0 = d_1 + d_3;                                  const double tmp0 = w12*(d_0 + d_1 + d_2 + d_3);
8322                                  const double tmp0_0 = d_0 + d_1 + d_2 + d_3;                                  const double tmp1 = w10*(d_1 + d_3);
8323                                  const double tmp1_0 = d_0 + d_2;                                  const double tmp2 = w11*(d_0 + d_2);
8324                                  const double tmp3_0 = d_0 + d_1;                                  const double tmp3 = w10*(d_2 + d_3);
8325                                  const double tmp6_0 = d_0 + d_3;                                  const double tmp4 = w11*(d_0 + d_1);
8326                                  const double tmp4_0 = d_2 + d_3;                                  const double tmp5 = w10*(d_0 + d_1);
8327                                  const double tmp5_0 = d_1 + d_2;                                  const double tmp6 = w11*(d_2 + d_3);
8328                                  const double tmp1_1 = tmp1_0*w25;                                  const double tmp7 = w12*(d_1 + d_2);
8329                                  const double tmp11_1 = d_0*w27;                                  const double tmp8 = w10*(d_0 + d_2);
8330                                  const double tmp9_1 = tmp5_0*w26;                                  const double tmp9 = w11*(d_1 + d_3);
8331                                  const double tmp12_1 = tmp2_0*w25;                                  const double tmp10 = w12*(d_0 + d_3);
8332                                  const double tmp3_1 = tmp3_0*w25;                                  EM_S[INDEX2(4,4,8)]+=d_0*w14 + d_3*w13 + tmp7;
8333                                  const double tmp6_1 = tmp3_0*w24;                                  EM_S[INDEX2(4,5,8)]+=tmp3 + tmp4;
8334                                  const double tmp2_1 = tmp2_0*w24;                                  EM_S[INDEX2(4,6,8)]+=tmp1 + tmp2;
8335                                  const double tmp10_1 = d_3*w28;                                  EM_S[INDEX2(4,7,8)]+=tmp0;
8336                                  const double tmp16_1 = tmp6_0*w26;                                  EM_S[INDEX2(5,4,8)]+=tmp3 + tmp4;
8337                                  const double tmp17_1 = d_1*w28;                                  EM_S[INDEX2(5,5,8)]+=d_1*w14 + d_2*w13 + tmp10;
8338                                  const double tmp7_1 = d_0*w28;                                  EM_S[INDEX2(5,6,8)]+=tmp0;
8339                                  const double tmp14_1 = d_2*w28;                                  EM_S[INDEX2(5,7,8)]+=tmp8 + tmp9;
8340                                  const double tmp18_1 = d_2*w27;                                  EM_S[INDEX2(6,4,8)]+=tmp1 + tmp2;
8341                                  const double tmp15_1 = d_1*w27;                                  EM_S[INDEX2(6,5,8)]+=tmp0;
8342                                  const double tmp0_1 = tmp0_0*w26;                                  EM_S[INDEX2(6,6,8)]+=d_1*w13 + d_2*w14 + tmp10;
8343                                  const double tmp4_1 = tmp4_0*w24;                                  EM_S[INDEX2(6,7,8)]+=tmp5 + tmp6;
8344                                  const double tmp5_1 = tmp4_0*w25;                                  EM_S[INDEX2(7,4,8)]+=tmp0;
8345                                  const double tmp8_1 = d_3*w27;                                  EM_S[INDEX2(7,5,8)]+=tmp8 + tmp9;
8346                                  const double tmp13_1 = tmp1_0*w24;                                  EM_S[INDEX2(7,6,8)]+=tmp5 + tmp6;
8347                                  EM_S[INDEX2(4,4,8)]+=tmp7_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX2(7,7,8)]+=d_0*w13 + d_3*w14 + tmp7;
8348                                  EM_S[INDEX2(5,4,8)]+=tmp3_1 + tmp4_1;                              } else { // constant data
8349                                  EM_S[INDEX2(6,4,8)]+=tmp1_1 + tmp2_1;                                  const double d_0 = d_p[0];
8350                                  EM_S[INDEX2(7,4,8)]+=tmp0_1;                                  EM_S[INDEX2(4,4,8)]+=16*d_0*w12;
8351                                  EM_S[INDEX2(4,5,8)]+=tmp3_1 + tmp4_1;                                  EM_S[INDEX2(4,5,8)]+=8*d_0*w12;
8352                                  EM_S[INDEX2(5,5,8)]+=tmp16_1 + tmp17_1 + tmp18_1;                                  EM_S[INDEX2(4,6,8)]+=8*d_0*w12;
8353                                  EM_S[INDEX2(6,5,8)]+=tmp0_1;                                  EM_S[INDEX2(4,7,8)]+=4*d_0*w12;
8354                                  EM_S[INDEX2(7,5,8)]+=tmp12_1 + tmp13_1;                                  EM_S[INDEX2(5,4,8)]+=8*d_0*w12;
8355                                  EM_S[INDEX2(4,6,8)]+=tmp1_1 + tmp2_1;                                  EM_S[INDEX2(5,5,8)]+=16*d_0*w12;
8356                                  EM_S[INDEX2(5,6,8)]+=tmp0_1;                                  EM_S[INDEX2(5,6,8)]+=4*d_0*w12;
8357                                  EM_S[INDEX2(6,6,8)]+=tmp14_1 + tmp15_1 + tmp16_1;                                  EM_S[INDEX2(5,7,8)]+=8*d_0*w12;
8358                                  EM_S[INDEX2(7,6,8)]+=tmp5_1 + tmp6_1;                                  EM_S[INDEX2(6,4,8)]+=8*d_0*w12;
8359                                  EM_S[INDEX2(4,7,8)]+=tmp0_1;                                  EM_S[INDEX2(6,5,8)]+=4*d_0*w12;
8360                                  EM_S[INDEX2(5,7,8)]+=tmp12_1 + tmp13_1;                                  EM_S[INDEX2(6,6,8)]+=16*d_0*w12;
8361                                  EM_S[INDEX2(6,7,8)]+=tmp5_1 + tmp6_1;                                  EM_S[INDEX2(6,7,8)]+=8*d_0*w12;
8362                                  EM_S[INDEX2(7,7,8)]+=tmp10_1 + tmp11_1 + tmp9_1;                                  EM_S[INDEX2(7,4,8)]+=4*d_0*w12;
8363                              } else { /* constant data */                                  EM_S[INDEX2(7,5,8)]+=8*d_0*w12;
8364                                  const double tmp2_1 = d_p[0]*w31;                                  EM_S[INDEX2(7,6,8)]+=8*d_0*w12;
8365                                  const double tmp0_1 = d_p[0]*w30;                                  EM_S[INDEX2(7,7,8)]+=16*d_0*w12;
                                 const double tmp1_1 = d_p[0]*w29;  
                                 EM_S[INDEX2(4,4,8)]+=tmp2_1;  
                                 EM_S[INDEX2(5,4,8)]+=tmp1_1;  
                                 EM_S[INDEX2(6,4,8)]+=tmp1_1;  
                                 EM_S[INDEX2(7,4,8)]+=tmp0_1;  
                                 EM_S[INDEX2(4,6,8)]+=tmp1_1;  
                                 EM_S[INDEX2(5,6,8)]+=tmp0_1;  
                                 EM_S[INDEX2(6,6,8)]+=tmp2_1;  
                                 EM_S[INDEX2(7,6,8)]+=tmp1_1;  
                                 EM_S[INDEX2(4,5,8)]+=tmp1_1;  
                                 EM_S[INDEX2(5,5,8)]+=tmp2_1;  
                                 EM_S[INDEX2(6,5,8)]+=tmp0_1;  
                                 EM_S[INDEX2(7,5,8)]+=tmp1_1;  
                                 EM_S[INDEX2(4,7,8)]+=tmp0_1;  
                                 EM_S[INDEX2(5,7,8)]+=tmp1_1;  
                                 EM_S[INDEX2(6,7,8)]+=tmp1_1;  
                                 EM_S[INDEX2(7,7,8)]+=tmp2_1;  
8366                              }                              }
8367                          }                          }
8368                          ///////////////                          ///////////////
# Line 8567  void Brick::assemblePDEBoundarySingle(Pa Line 8375  void Brick::assemblePDEBoundarySingle(Pa
8375                                  const double y_1 = y_p[1];                                  const double y_1 = y_p[1];
8376                                  const double y_2 = y_p[2];                                  const double y_2 = y_p[2];
8377                                  const double y_3 = y_p[3];                                  const double y_3 = y_p[3];
8378                                  const double tmp0_0 = y_1 + y_2;                                  const double tmp0 = 6*w12*(y_1 + y_2);
8379                                  const double tmp1_0 = y_0 + y_3;                                  const double tmp1 = 6*w12*(y_0 + y_3);
8380                                  const double tmp7_1 = w34*y_1;                                  EM_F[4]+=tmp0 + 6*w10*y_3 + 6*w11*y_0;
8381                                  const double tmp3_1 = w32*y_1;                                  EM_F[5]+=tmp1 + 6*w10*y_2 + 6*w11*y_1;
8382                                  const double tmp8_1 = w32*y_3;                                  EM_F[6]+=tmp1 + 6*w10*y_1 + 6*w11*y_2;
8383                                  const double tmp0_1 = w32*y_0;                                  EM_F[7]+=tmp0 + 6*w10*y_0 + 6*w11*y_3;
8384                                  const double tmp2_1 = w34*y_3;                              } else { // constant data
8385                                  const double tmp9_1 = w34*y_0;                                  EM_F[4]+=36*w12*y_p[0];
8386                                  const double tmp6_1 = w32*y_2;                                  EM_F[5]+=36*w12*y_p[0];
8387                                  const double tmp5_1 = w34*y_2;                                  EM_F[6]+=36*w12*y_p[0];
8388                                  const double tmp1_1 = tmp0_0*w33;                                  EM_F[7]+=36*w12*y_p[0];
                                 const double tmp4_1 = tmp1_0*w33;  
                                 EM_F[4]+=tmp0_1 + tmp1_1 + tmp2_1;  
                                 EM_F[5]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                 EM_F[6]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                 EM_F[7]+=tmp1_1 + tmp8_1 + tmp9_1;  
                             } else { /* constant data */  
                                 const double tmp0_1 = w35*y_p[0];  
                                 EM_F[4]+=tmp0_1;  
                                 EM_F[5]+=tmp0_1;  
                                 EM_F[6]+=tmp0_1;  
                                 EM_F[7]+=tmp0_1;  
8389                              }                              }
8390                          }                          }
8391                          const index_t firstNode=m_NN[0]*m_NN[1]*(m_NN[2]-2)+m_NN[0]*k1+k0;                          const index_t firstNode=m_NN[0]*m_NN[1]*(m_NN[2]-2)+m_NN[0]*k1+k0;
# Line 8605  void Brick::assemblePDEBoundarySingle(Pa Line 8402  void Brick::assemblePDEBoundarySingle(Pa
8402  void Brick::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,  void Brick::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
8403        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
8404  {  {
8405      const double w0 = 0.0625*m_dx[0]*m_dx[1];      const double w0 = m_dx[0]*m_dx[1]/16.;
8406      const double w1 = 0.0625*m_dx[0]*m_dx[2];      const double w1 = m_dx[0]*m_dx[2]/16.;
8407      const double w2 = 0.0625*m_dx[1]*m_dx[2];      const double w2 = m_dx[1]*m_dx[2]/16.;
8408      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
8409      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
8410      rhs.requireWrite();      rhs.requireWrite();
# Line 8907  void Brick::assemblePDEBoundarySingleRed Line 8704  void Brick::assemblePDEBoundarySingleRed
8704  void Brick::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,  void Brick::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
8705        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
8706  {  {
     const double h0 = m_dx[0];  
     const double h1 = m_dx[1];  
     const double h2 = m_dx[2];  
8707      dim_t numEq, numComp;      dim_t numEq, numComp;
8708      if (!mat)      if (!mat)
8709          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
# Line 8917  void Brick::assemblePDEBoundarySystem(Pa Line 8711  void Brick::assemblePDEBoundarySystem(Pa
8711          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
8712          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
8713      }      }
8714      const double w0 = 0.0018607582807716854616*h1*h2;      /* GENERATOR SNIP_PDEBC_SYSTEM_PRE TOP */
8715      const double w1 = 0.025917019497006092316*h1*h2;      const double SQRT3 = 1.73205080756887719318;
8716      const double w10 = 0.01116454968463011277*h1*h2;      const double w12 = m_dx[0]*m_dx[1]/144;
8717      const double w11 = 0.25*h1*h2;      const double w11 = w12*(SQRT3 + 2);
8718      const double w12 = 0.0018607582807716854616*h0*h2;      const double w14 = w12*(4*SQRT3 + 7);
8719      const double w13 = 0.025917019497006092316*h0*h2;      const double w13 = w12*(-4*SQRT3 + 7);
8720      const double w14 = 0.0069444444444444444444*h0*h2;      const double w10 = w12*(-SQRT3 + 2);
8721      const double w15 = 0.00049858867864229740201*h0*h2;      const double w7 = m_dx[0]*m_dx[2]/144;
8722      const double w16 = 0.09672363354357992482*h0*h2;      const double w6 = w7*(SQRT3 + 2);
8723      const double w17 = 0.055555555555555555556*h0*h2;      const double w8 = w7*(-4*SQRT3 + 7);
8724      const double w18 = 0.11111111111111111111*h0*h2;      const double w9 = w7*(4*SQRT3 + 7);
8725      const double w19 = 0.027777777777777777778*h0*h2;      const double w5 = w7*(-SQRT3 + 2);
8726      const double w2 = 0.0069444444444444444444*h1*h2;      const double w2 = m_dx[1]*m_dx[2]/144;
8727      const double w20 = 0.1555021169820365539*h0*h2;      const double w1 = w2*(SQRT3 + 2);
8728      const double w21 = 0.041666666666666666667*h0*h2;      const double w4 = w2*(4*SQRT3 + 7);
8729      const double w22 = 0.01116454968463011277*h0*h2;      const double w3 = w2*(-4*SQRT3 + 7);
8730      const double w23 = 0.25*h0*h2;      const double w0 = w2*(-SQRT3 + 2);
8731      const double w24 = 0.0018607582807716854616*h0*h1;      /* GENERATOR SNIP_PDEBC_SYSTEM_PRE BOTTOM */
     const double w25 = 0.025917019497006092316*h0*h1;  
     const double w26 = 0.0069444444444444444444*h0*h1;  
     const double w27 = 0.00049858867864229740201*h0*h1;  
     const double w28 = 0.09672363354357992482*h0*h1;  
     const double w29 = 0.055555555555555555556*h0*h1;  
     const double w3 = 0.00049858867864229740201*h1*h2;  
     const double w30 = 0.027777777777777777778*h0*h1;  
     const double w31 = 0.11111111111111111111*h0*h1;  
     const double w32 = 0.1555021169820365539*h0*h1;  
     const double w33 = 0.041666666666666666667*h0*h1;  
     const double w34 = 0.01116454968463011277*h0*h1;  
     const double w35 = 0.25*h0*h1;  
     const double w4 = 0.09672363354357992482*h1*h2;  
     const double w5 = 0.055555555555555555556*h1*h2;  
     const double w6 = 0.11111111111111111111*h1*h2;  
     const double w7 = 0.027777777777777777778*h1*h2;  
     const double w8 = 0.1555021169820365539*h1*h2;  
     const double w9 = 0.041666666666666666667*h1*h2;  
8732      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
8733      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
8734      rhs.requireWrite();      rhs.requireWrite();
# Line 8966  void Brick::assemblePDEBoundarySystem(Pa Line 8742  void Brick::assemblePDEBoundarySystem(Pa
8742                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
8743                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
8744                          const index_t e = INDEX2(k1,k2,m_NE[1]);                          const index_t e = INDEX2(k1,k2,m_NE[1]);
8745                            /* GENERATOR SNIP_PDEBC_SYSTEM_0 TOP */
8746                          ///////////////                          ///////////////
8747                          // process d //                          // process d //
8748                          ///////////////                          ///////////////
# Line 8974  void Brick::assemblePDEBoundarySystem(Pa Line 8751  void Brick::assemblePDEBoundarySystem(Pa
8751                              if (d.actsExpanded()) {                              if (d.actsExpanded()) {
8752                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8753                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8754                                          const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                          const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
8755                                          const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                          const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
8756                                          const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];                                          const double d_2 = d_p[INDEX3(k,m,2,numEq,numComp)];
8757                                          const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];                                          const double d_3 = d_p[INDEX3(k,m,3,numEq,numComp)];
8758                                          const double tmp3_0 = d_1 + d_3;                                          const double tmp0 = w0*(d_0 + d_1);
8759                                          const double tmp6_0 = d_0 + d_1 + d_2 + d_3;                                          const double tmp1 = w1*(d_2 + d_3);
8760                                          const double tmp2_0 = d_0 + d_2;                                          const double tmp2 = w0*(d_0 + d_2);
8761                                          const double tmp0_0 = d_0 + d_1;                                          const double tmp3 = w1*(d_1 + d_3);
8762                                          const double tmp4_0 = d_0 + d_3;                                          const double tmp4 = w0*(d_1 + d_3);
8763                                          const double tmp1_0 = d_2 + d_3;                                          const double tmp5 = w1*(d_0 + d_2);
8764                                          const double tmp5_0 = d_1 + d_2;                                          const double tmp6 = w0*(d_2 + d_3);
8765                                          const double tmp14_1 = d_0*w4;                                          const double tmp7 = w1*(d_0 + d_1);
8766                                          const double tmp17_1 = d_3*w4;                                          const double tmp8 = w2*(d_0 + d_3);
8767                                          const double tmp12_1 = d_1*w4;                                          const double tmp9 = w2*(d_1 + d_2);
8768                                          const double tmp9_1 = d_2*w4;                                          const double tmp10 = w2*(d_0 + d_1 + d_2 + d_3);
8769                                          const double tmp4_1 = tmp3_0*w0;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=d_0*w4 + d_3*w3 + tmp9;
8770                                          const double tmp3_1 = tmp3_0*w1;                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp6 + tmp7;
8771                                          const double tmp7_1 = tmp0_0*w1;                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp4 + tmp5;
8772                                          const double tmp8_1 = d_1*w3;                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp10;
8773                                          const double tmp16_1 = d_0*w3;                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp6 + tmp7;
8774                                          const double tmp18_1 = tmp6_0*w2;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=d_1*w4 + d_2*w3 + tmp8;
8775                                          const double tmp1_1 = tmp1_0*w1;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp10;
8776                                          const double tmp15_1 = tmp5_0*w2;                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp2 + tmp3;
8777                                          const double tmp6_1 = tmp1_0*w0;                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp4 + tmp5;
8778                                          const double tmp2_1 = tmp2_0*w0;                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp10;
8779                                          const double tmp13_1 = d_3*w3;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=d_1*w3 + d_2*w4 + tmp8;
8780                                          const double tmp5_1 = tmp2_0*w1;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0 + tmp1;
8781                                          const double tmp10_1 = tmp4_0*w2;                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp10;
8782                                          const double tmp11_1 = d_2*w3;                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp2 + tmp3;
8783                                          const double tmp0_1 = tmp0_0*w0;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0 + tmp1;
8784                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0_1 + tmp1_1;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=d_0*w3 + d_3*w4 + tmp9;
                                         EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp2_1 + tmp3_1;  
                                         EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp4_1 + tmp5_1;  
                                         EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp6_1 + tmp7_1;  
                                         EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp10_1 + tmp8_1 + tmp9_1;  
                                         EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp12_1;  
                                         EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp13_1 + tmp14_1 + tmp15_1;  
                                         EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp15_1 + tmp16_1 + tmp17_1;  
                                         EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp4_1 + tmp5_1;  
                                         EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp18_1;  
                                         EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp18_1;  
                                         EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1 + tmp1_1;  
                                         EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp6_1 + tmp7_1;  
                                         EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp18_1;  
                                         EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp2_1 + tmp3_1;  
                                         EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp18_1;  
8785                                      }                                      }
8786                                   }                                   }
8787                              } else { /* constant data */                              } else { // constant data
8788                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8789                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8790                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double d_0 = d_p[INDEX2(k, m, numEq)];
8791                                          const double tmp0_1 = d_0*w5;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=16*d_0*w2;
8792                                          const double tmp2_1 = d_0*w7;                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=8*d_0*w2;
8793                                          const double tmp1_1 = d_0*w6;                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=8*d_0*w2;
8794                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=4*d_0*w2;
8795                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=8*d_0*w2;
8796                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=16*d_0*w2;
8797                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=4*d_0*w2;
8798                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=8*d_0*w2;
8799                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=8*d_0*w2;
8800                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=4*d_0*w2;
8801                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=16*d_0*w2;
8802                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=8*d_0*w2;
8803                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=4*d_0*w2;
8804                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=8*d_0*w2;
8805                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=8*d_0*w2;
8806                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=16*d_0*w2;
                                         EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp2_1;  
                                         EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp2_1;  
8807                                      }                                      }
8808                                  }                                  }
8809                              }                              }
# Line 9060  void Brick::assemblePDEBoundarySystem(Pa Line 8819  void Brick::assemblePDEBoundarySystem(Pa
8819                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];
8820                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];
8821                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];
8822                                      const double tmp0_0 = y_1 + y_2;                                      const double tmp0 = 6*w2*(y_1 + y_2);
8823                                      const double tmp1_0 = y_0 + y_3;                                      const double tmp1 = 6*w2*(y_0 + y_3);
8824                                      const double tmp2_1 = w10*y_3;                                      EM_F[INDEX2(k,0,numEq)]+=tmp0 + 6*w0*y_3 + 6*w1*y_0;
8825                                      const double tmp8_1 = w10*y_0;                                      EM_F[INDEX2(k,2,numEq)]+=tmp1 + 6*w0*y_2 + 6*w1*y_1;
8826                                      const double tmp5_1 = w10*y_2;                                      EM_F[INDEX2(k,4,numEq)]+=tmp1 + 6*w0*y_1 + 6*w1*y_2;
8827                                      const double tmp3_1 = w8*y_1;                                      EM_F[INDEX2(k,6,numEq)]+=tmp0 + 6*w0*y_0 + 6*w1*y_3;
8828                                      const double tmp9_1 = w8*y_3;                                  }
8829                                      const double tmp0_1 = w8*y_0;                              } else { // constant data
8830                                      const double tmp1_1 = tmp0_0*w9;                                  for (index_t k=0; k<numEq; k++) {
8831                                      const double tmp7_1 = w8*y_2;                                      EM_F[INDEX2(k,0,numEq)]+=36*w2*y_p[k];
8832                                      const double tmp4_1 = tmp1_0*w9;                                      EM_F[INDEX2(k,2,numEq)]+=36*w2*y_p[k];
8833                                      const double tmp6_1 = w10*y_1;                                      EM_F[INDEX2(k,4,numEq)]+=36*w2*y_p[k];
8834                                      EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_F[INDEX2(k,6,numEq)]+=36*w2*y_p[k];
                                     EM_F[INDEX2(k,2,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_F[INDEX2(k,4,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                     EM_F[INDEX2(k,6,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;  
                                 }  
                             } else { /* constant data */  
                                 for (index_t k=0; k<numEq; k++) {  
                                     const double y_0 = y_p[k];  
                                     const double tmp0_1 = w11*y_0;  
                                     EM_F[INDEX2(k,0,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,2,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,4,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,6,numEq)]+=tmp0_1;  
8835                                  }                                  }
8836                              }                              }
8837                          }                          }
8838                            /* GENERATOR SNIP_PDEBC_SYSTEM_0 BOTTOM */
8839                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*k1;                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*k1;
8840                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
8841                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 9104  void Brick::assemblePDEBoundarySystem(Pa Line 8852  void Brick::assemblePDEBoundarySystem(Pa
8852                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
8853                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
8854                          const index_t e = m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]);                          const index_t e = m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]);
8855                            /* GENERATOR SNIP_PDEBC_SYSTEM_1 TOP */
8856                          ///////////////                          ///////////////
8857                          // process d //                          // process d //
8858                          ///////////////                          ///////////////
# Line 9112  void Brick::assemblePDEBoundarySystem(Pa Line 8861  void Brick::assemblePDEBoundarySystem(Pa
8861                              if (d.actsExpanded()) {                              if (d.actsExpanded()) {
8862                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8863                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8864                                          const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                          const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
8865                                          const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                          const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
8866                                          const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];                                          const double d_2 = d_p[INDEX3(k,m,2,numEq,numComp)];
8867                                          const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];                                          const double d_3 = d_p[INDEX3(k,m,3,numEq,numComp)];
8868                                          const double tmp1_0 = d_1 + d_3;                                          const double tmp0 = w0*(d_0 + d_2);
8869                                          const double tmp6_0 = d_0 + d_1 + d_2 + d_3;                                          const double tmp1 = w1*(d_1 + d_3);
8870                                          const double tmp0_0 = d_0 + d_2;                                          const double tmp2 = w0*(d_2 + d_3);
8871                                          const double tmp3_0 = d_0 + d_1;                                          const double tmp3 = w1*(d_0 + d_1);
8872                                          const double tmp4_0 = d_0 + d_3;                                          const double tmp4 = w0*(d_1 + d_3);
8873                                          const double tmp2_0 = d_2 + d_3;                                          const double tmp5 = w1*(d_0 + d_2);
8874                                          const double tmp5_0 = d_1 + d_2;                                          const double tmp6 = w2*(d_0 + d_3);
8875                                          const double tmp10_1 = d_3*w4;                                          const double tmp7 = w2*(d_1 + d_2);
8876                                          const double tmp13_1 = tmp2_0*w1;                                          const double tmp8 = w0*(d_0 + d_1);
8877                                          const double tmp16_1 = d_0*w4;                                          const double tmp9 = w1*(d_2 + d_3);
8878                                          const double tmp18_1 = d_2*w4;                                          const double tmp10 = w2*(d_0 + d_1 + d_2 + d_3);
8879                                          const double tmp12_1 = tmp3_0*w0;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=d_0*w4 + d_3*w3 + tmp7;
8880                                          const double tmp3_1 = tmp3_0*w1;                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp2 + tmp3;
8881                                          const double tmp5_1 = tmp0_0*w1;                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp4 + tmp5;
8882                                          const double tmp17_1 = d_1*w3;                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp10;
8883                                          const double tmp9_1 = d_0*w3;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp2 + tmp3;
8884                                          const double tmp14_1 = tmp6_0*w2;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=d_1*w4 + d_2*w3 + tmp6;
8885                                          const double tmp1_1 = tmp1_0*w1;                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp10;
8886                                          const double tmp11_1 = tmp5_0*w2;                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0 + tmp1;
8887                                          const double tmp4_1 = tmp1_0*w0;                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp4 + tmp5;
8888                                          const double tmp2_1 = tmp2_0*w0;                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp10;
8889                                          const double tmp15_1 = d_3*w3;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=d_1*w3 + d_2*w4 + tmp6;
8890                                          const double tmp7_1 = d_1*w4;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp8 + tmp9;
8891                                          const double tmp8_1 = tmp4_0*w2;                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp10;
8892                                          const double tmp6_1 = d_2*w3;                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0 + tmp1;
8893                                          const double tmp0_1 = tmp0_0*w0;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp8 + tmp9;
8894                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=d_0*w3 + d_3*w4 + tmp7;
                                         EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp2_1 + tmp3_1;  
                                         EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp4_1 + tmp5_1;  
                                         EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1 + tmp1_1;  
                                         EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp6_1 + tmp7_1 + tmp8_1;  
                                         EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp4_1 + tmp5_1;  
                                         EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp9_1;  
                                         EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp12_1 + tmp13_1;  
                                         EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp14_1;  
                                         EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp11_1 + tmp15_1 + tmp16_1;  
                                         EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp14_1;  
                                         EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp17_1 + tmp18_1 + tmp8_1;  
                                         EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp12_1 + tmp13_1;  
                                         EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp14_1;  
                                         EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp2_1 + tmp3_1;  
                                         EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp14_1;  
8895                                      }                                      }
8896                                   }                                   }
8897                              } else { /* constant data */                              } else { // constant data
8898                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8899                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8900                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double d_0 = d_p[INDEX2(k, m, numEq)];
8901                                          const double tmp0_1 = d_0*w5;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=16*d_0*w2;
8902                                          const double tmp2_1 = d_0*w7;                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=8*d_0*w2;
8903                                          const double tmp1_1 = d_0*w6;                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=8*d_0*w2;
8904                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=4*d_0*w2;
8905                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=8*d_0*w2;
8906                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=16*d_0*w2;
8907                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=4*d_0*w2;
8908                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=8*d_0*w2;
8909                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=8*d_0*w2;
8910                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=4*d_0*w2;
8911                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=16*d_0*w2;
8912                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=8*d_0*w2;
8913                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=4*d_0*w2;
8914                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=8*d_0*w2;
8915                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=8*d_0*w2;
8916                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=16*d_0*w2;
                                         EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp2_1;  
                                         EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp2_1;  
8917                                      }                                      }
8918                                  }                                  }
8919                              }                              }
# Line 9198  void Brick::assemblePDEBoundarySystem(Pa Line 8929  void Brick::assemblePDEBoundarySystem(Pa
8929                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];
8930                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];
8931                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];
8932                                      const double tmp0_0 = y_1 + y_2;                                      const double tmp0 = 6*w2*(y_1 + y_2);
8933                                      const double tmp1_0 = y_0 + y_3;                                      const double tmp1 = 6*w2*(y_0 + y_3);
8934                                      const double tmp2_1 = w10*y_3;                                      EM_F[INDEX2(k,1,numEq)]+=tmp0 + 6*w0*y_3 + 6*w1*y_0;
8935                                      const double tmp8_1 = w10*y_0;                                      EM_F[INDEX2(k,3,numEq)]+=tmp1 + 6*w0*y_2 + 6*w1*y_1;
8936                                      const double tmp5_1 = w10*y_2;                                      EM_F[INDEX2(k,5,numEq)]+=tmp1 + 6*w0*y_1 + 6*w1*y_2;
8937                                      const double tmp3_1 = w8*y_1;                                      EM_F[INDEX2(k,7,numEq)]+=tmp0 + 6*w0*y_0 + 6*w1*y_3;
8938                                      const double tmp9_1 = w8*y_3;                                  }
8939                                      const double tmp0_1 = w8*y_0;                              } else { // constant data
8940                                      const double tmp1_1 = tmp0_0*w9;                                  for (index_t k=0; k<numEq; k++) {
8941                                      const double tmp7_1 = w8*y_2;                                      EM_F[INDEX2(k,1,numEq)]+=36*w2*y_p[k];
8942                                      const double tmp4_1 = tmp1_0*w9;                                      EM_F[INDEX2(k,3,numEq)]+=36*w2*y_p[k];
8943                                      const double tmp6_1 = w10*y_1;                                      EM_F[INDEX2(k,5,numEq)]+=36*w2*y_p[k];
8944                                      EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_F[INDEX2(k,7,numEq)]+=36*w2*y_p[k];
                                     EM_F[INDEX2(k,3,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_F[INDEX2(k,5,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                     EM_F[INDEX2(k,7,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;  
                                 }  
                             } else { /* constant data */  
                                 for (index_t k=0; k<numEq; k++) {  
                                     const double y_0 = y_p[k];  
                                     const double tmp0_1 = w11*y_0;  
                                     EM_F[INDEX2(k,1,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,5,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,7,numEq)]+=tmp0_1;  
8945                                  }                                  }
8946                              }                              }
8947                          }                          }
8948                            /* GENERATOR SNIP_PDEBC_SYSTEM_1 BOTTOM */
8949                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*(k1+1)-2;                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*(k1+1)-2;
8950                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
8951                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 9242  void Brick::assemblePDEBoundarySystem(Pa Line 8962  void Brick::assemblePDEBoundarySystem(Pa
8962                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
8963                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
8964                          const index_t e = m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]);                          const index_t e = m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]);
8965                            /* GENERATOR SNIP_PDEBC_SYSTEM_2 TOP */
8966                          ///////////////                          ///////////////
8967                          // process d //                          // process d //
8968                          ///////////////                          ///////////////
# Line 9250  void Brick::assemblePDEBoundarySystem(Pa Line 8971  void Brick::assemblePDEBoundarySystem(Pa
8971                              if (d.actsExpanded()) {                              if (d.actsExpanded()) {
8972                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8973                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8974                                          const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                          const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
8975                                          const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                          const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
8976                                          const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];                                          const double d_2 = d_p[INDEX3(k,m,2,numEq,numComp)];
8977                                          const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];                                          const double d_3 = d_p[INDEX3(k,m,3,numEq,numComp)];
8978                                          const double tmp2_0 = d_1 + d_3;                                          const double tmp0 = w5*(d_0 + d_1);
8979                                          const double tmp5_0 = d_0 + d_1 + d_2 + d_3;                                          const double tmp1 = w6*(d_2 + d_3);
8980                                          const double tmp3_0 = d_0 + d_2;                                          const double tmp2 = w5*(d_0 + d_2);
8981                                          const double tmp1_0 = d_0 + d_1;                                          const double tmp3 = w6*(d_1 + d_3);
8982                                          const double tmp4_0 = d_0 + d_3;                                          const double tmp4 = w5*(d_1 + d_3);
8983                                          const double tmp0_0 = d_2 + d_3;                                          const double tmp5 = w6*(d_0 + d_2);
8984                                          const double tmp6_0 = d_1 + d_2;                                          const double tmp6 = w7*(d_0 + d_3);
8985                                          const double tmp2_1 = tmp2_0*w13;                                          const double tmp7 = w7*(d_0 + d_1 + d_2 + d_3);
8986                                          const double tmp14_1 = d_3*w15;                                          const double tmp8 = w7*(d_1 + d_2);
8987                                          const double tmp0_1 = tmp0_0*w13;                                          const double tmp9 = w5*(d_2 + d_3);
8988                                          const double tmp3_1 = tmp3_0*w12;                                          const double tmp10 = w6*(d_0 + d_1);
8989                                          const double tmp17_1 = tmp1_0*w13;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=d_0*w9 + d_3*w8 + tmp8;
8990                                          const double tmp18_1 = tmp0_0*w12;                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp10 + tmp9;
8991                                          const double tmp8_1 = d_1*w15;                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp4 + tmp5;
8992                                          const double tmp16_1 = d_0*w15;                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp7;
8993                                          const double tmp11_1 = d_2*w15;                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp10 + tmp9;
8994                                          const double tmp5_1 = tmp2_0*w12;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=d_1*w9 + d_2*w8 + tmp6;
8995                                          const double tmp15_1 = d_3*w16;                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp7;
8996                                          const double tmp13_1 = tmp6_0*w14;                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp2 + tmp3;
8997                                          const double tmp1_1 = tmp1_0*w12;                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp4 + tmp5;
8998                                          const double tmp7_1 = tmp4_0*w14;                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp7;
8999                                          const double tmp12_1 = d_0*w16;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=d_1*w8 + d_2*w9 + tmp6;
9000                                          const double tmp10_1 = d_1*w16;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0 + tmp1;
9001                                          const double tmp6_1 = d_2*w16;                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp7;
9002                                          const double tmp9_1 = tmp5_0*w14;                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp2 + tmp3;
9003                                          const double tmp4_1 = tmp3_0*w13;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0 + tmp1;
9004                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0_1 + tmp1_1;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=d_0*w8 + d_3*w9 + tmp8;
                                         EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp2_1 + tmp3_1;  
                                         EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp4_1 + tmp5_1;  
                                         EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp6_1 + tmp7_1 + tmp8_1;  
                                         EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1;  
                                         EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp9_1;  
                                         EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp7_1;  
                                         EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp12_1 + tmp13_1 + tmp14_1;  
                                         EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp9_1;  
                                         EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0_1 + tmp1_1;  
                                         EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp4_1 + tmp5_1;  
                                         EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp13_1 + tmp15_1 + tmp16_1;  
                                         EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp9_1;  
                                         EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp17_1 + tmp18_1;  
                                         EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp17_1 + tmp18_1;  
                                         EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp9_1;  
9005                                      }                                      }
9006                                   }                                   }
9007                              } else { /* constant data */                              } else { // constant data
9008                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9009                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9010                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double d_0 = d_p[INDEX2(k, m, numEq)];
9011                                          const double tmp0_1 = d_0*w17;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=16*d_0*w7;
9012                                          const double tmp2_1 = d_0*w19;                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=8*d_0*w7;
9013                                          const double tmp1_1 = d_0*w18;                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=8*d_0*w7;
9014                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=4*d_0*w7;
9015                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=8*d_0*w7;
9016                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=16*d_0*w7;
9017                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=4*d_0*w7;
9018                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=8*d_0*w7;
9019                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=8*d_0*w7;
9020                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=4*d_0*w7;
9021                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=16*d_0*w7;
9022                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=8*d_0*w7;
9023                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=4*d_0*w7;
9024                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=8*d_0*w7;
9025                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=8*d_0*w7;
9026                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=16*d_0*w7;
                                         EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp2_1;  
9027                                      }                                      }
9028                                  }                                  }
9029                              }                              }
# Line 9336  void Brick::assemblePDEBoundarySystem(Pa Line 9039  void Brick::assemblePDEBoundarySystem(Pa
9039                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];
9040                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];
9041                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];
9042                                      const double tmp0_0 = y_1 + y_2;                                      const double tmp0 = 6*w7*(y_1 + y_2);
9043                                      const double tmp1_0 = y_0 + y_3;                                      const double tmp1 = 6*w7*(y_0 + y_3);
9044                                      const double tmp0_1 = w22*y_3;                                      EM_F[INDEX2(k,0,numEq)]+=tmp0 + 6*w5*y_3 + 6*w6*y_0;
9045                                      const double tmp6_1 = w22*y_1;                                      EM_F[INDEX2(k,1,numEq)]+=tmp1 + 6*w5*y_2 + 6*w6*y_1;
9046                                      const double tmp3_1 = w22*y_2;                                      EM_F[INDEX2(k,4,numEq)]+=tmp1 + 6*w5*y_1 + 6*w6*y_2;
9047                                      const double tmp5_1 = w20*y_1;                                      EM_F[INDEX2(k,5,numEq)]+=tmp0 + 6*w5*y_0 + 6*w6*y_3;
9048                                      const double tmp9_1 = w20*y_3;                                  }
9049                                      const double tmp4_1 = tmp1_0*w21;                              } else { // constant data
9050                                      const double tmp8_1 = w22*y_0;                                  for (index_t k=0; k<numEq; k++) {
9051                                      const double tmp2_1 = w20*y_0;                                      EM_F[INDEX2(k,0,numEq)]+=36*w7*y_p[k];
9052                                      const double tmp7_1 = w20*y_2;                                      EM_F[INDEX2(k,1,numEq)]+=36*w7*y_p[k];
9053                                      const double tmp1_1 = tmp0_0*w21;                                      EM_F[INDEX2(k,4,numEq)]+=36*w7*y_p[k];
9054                                      EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_F[INDEX2(k,5,numEq)]+=36*w7*y_p[k];
                                     EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_F[INDEX2(k,4,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                     EM_F[INDEX2(k,5,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;  
                                 }  
                             } else { /* constant data */  
                                 for (index_t k=0; k<numEq; k++) {  
                                     const double y_0 = y_p[k];  
                                     const double tmp0_1 = w23*y_0;  
                                     EM_F[INDEX2(k,0,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,1,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,4,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,5,numEq)]+=tmp0_1;  
9055                                  }                                  }
9056                              }                              }
9057                          }                          }
9058                            /* GENERATOR SNIP_PDEBC_SYSTEM_2 BOTTOM */
9059                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+k0;                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+k0;
9060                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9061                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 9380  void Brick::assemblePDEBoundarySystem(Pa Line 9072  void Brick::assemblePDEBoundarySystem(Pa
9072                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
9073                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
9074                          const index_t e = m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]);                          const index_t e = m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]);
9075                            /* GENERATOR SNIP_PDEBC_SYSTEM_3 TOP */
9076                          ///////////////                          ///////////////
9077                          // process d //                          // process d //
9078                          ///////////////                          ///////////////
# Line 9388  void Brick::assemblePDEBoundarySystem(Pa Line 9081  void Brick::assemblePDEBoundarySystem(Pa
9081                              if (d.actsExpanded()) {                              if (d.actsExpanded()) {
9082                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9083                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9084                                          const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                          const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
9085                                          const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                          const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
9086                                          const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];                                          const double d_2 = d_p[INDEX3(k,m,2,numEq,numComp)];
9087                                          const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];                                          const double d_3 = d_p[INDEX3(k,m,3,numEq,numComp)];
9088                                          const double tmp0_0 = d_1 + d_3;                                          const double tmp0 = w5*(d_0 + d_2);
9089                                          const double tmp2_0 = d_0 + d_1 + d_2 + d_3;                                          const double tmp1 = w6*(d_1 + d_3);
9090                                          const double tmp1_0 = d_0 + d_2;                                          const double tmp2 = w5*(d_1 + d_3);
9091                                          const double tmp4_0 = d_0 + d_1;                                          const double tmp3 = w6*(d_0 + d_2);
9092                                          const double tmp5_0 = d_0 + d_3;                                          const double tmp4 = w7*(d_0 + d_1 + d_2 + d_3);
9093                                          const double tmp3_0 = d_2 + d_3;                                          const double tmp5 = w5*(d_0 + d_1);
9094                                          const double tmp6_0 = d_1 + d_2;                                          const double tmp6 = w6*(d_2 + d_3);
9095                                          const double tmp15_1 = tmp4_0*w13;                                          const double tmp7 = w7*(d_0 + d_3);
9096                                          const double tmp10_1 = d_0*w16;                                          const double tmp8 = w7*(d_1 + d_2);
9097                                          const double tmp6_1 = tmp4_0*w12;                                          const double tmp9 = w5*(d_2 + d_3);
9098                                          const double tmp16_1 = tmp3_0*w12;                                          const double tmp10 = w6*(d_0 + d_1);
9099                                          const double tmp0_1 = tmp0_0*w13;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=d_0*w9 + d_3*w8 + tmp8;
9100                                          const double tmp2_1 = tmp1_0*w13;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp10 + tmp9;
9101                                          const double tmp18_1 = d_1*w15;                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp2 + tmp3;
9102                                          const double tmp14_1 = d_0*w15;                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp4;
9103                                          const double tmp9_1 = d_2*w15;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp10 + tmp9;
9104                                          const double tmp4_1 = tmp2_0*w14;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=d_1*w9 + d_2*w8 + tmp7;
9105                                          const double tmp13_1 = d_3*w16;                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp4;
9106                                          const double tmp11_1 = tmp6_0*w14;                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0 + tmp1;
9107                                          const double tmp1_1 = tmp1_0*w12;                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp2 + tmp3;
9108                                          const double tmp12_1 = d_3*w15;                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp4;
9109                                          const double tmp3_1 = tmp0_0*w12;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=d_1*w8 + d_2*w9 + tmp7;
9110                                          const double tmp7_1 = d_1*w16;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp5 + tmp6;
9111                                          const double tmp17_1 = d_2*w16;                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp4;
9112                                          const double tmp8_1 = tmp5_0*w14;                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0 + tmp1;
9113                                          const double tmp5_1 = tmp3_0*w13;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp5 + tmp6;
9114                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=d_0*w8 + d_3*w9 + tmp8;
                                         EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp2_1 + tmp3_1;  
                                         EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1 + tmp1_1;  
                                         EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp4_1;  
                                         EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp5_1 + tmp6_1;  
                                         EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                         EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp5_1 + tmp6_1;  
                                         EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp4_1;  
                                         EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp4_1;  
                                         EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp12_1;  
                                         EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp11_1 + tmp13_1 + tmp14_1;  
                                         EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp4_1;  
                                         EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp15_1 + tmp16_1;  
                                         EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp17_1 + tmp18_1 + tmp8_1;  
                                         EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp15_1 + tmp16_1;  
                                         EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp2_1 + tmp3_1;  
9115                                      }                                      }
9116                                   }                                   }
9117                              } else { /* constant data */                              } else { // constant data
9118                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9119                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9120                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double d_0 = d_p[INDEX2(k, m, numEq)];
9121                                          const double tmp0_1 = d_0*w17;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=16*d_0*w7;
9122                                          const double tmp1_1 = d_0*w19;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=8*d_0*w7;
9123                                          const double tmp2_1 = d_0*w18;                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=8*d_0*w7;
9124                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=4*d_0*w7;
9125                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=8*d_0*w7;
9126                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=16*d_0*w7;
9127                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=4*d_0*w7;
9128                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=8*d_0*w7;
9129                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=8*d_0*w7;
9130                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=4*d_0*w7;
9131                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=16*d_0*w7;
9132                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=8*d_0*w7;
9133                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=4*d_0*w7;
9134                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=8*d_0*w7;
9135                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=8*d_0*w7;
9136                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=16*d_0*w7;
                                         EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp2_1;  
                                         EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1;  
9137                                      }                                      }
9138                                  }                                  }
9139                              }                              }
# Line 9474  void Brick::assemblePDEBoundarySystem(Pa Line 9149  void Brick::assemblePDEBoundarySystem(Pa
9149                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];
9150                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];
9151                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];
9152                                      const double tmp0_0 = y_1 + y_2;                                      const double tmp0 = 6*w7*(y_1 + y_2);
9153                                      const double tmp1_0 = y_0 + y_3;                                      const double tmp1 = 6*w7*(y_0 + y_3);
9154                                      const double tmp0_1 = w22*y_3;                                      EM_F[INDEX2(k,2,numEq)]+=tmp0 + 6*w5*y_3 + 6*w6*y_0;
9155                                      const double tmp6_1 = w22*y_1;                                      EM_F[INDEX2(k,3,numEq)]+=tmp1 + 6*w5*y_2 + 6*w6*y_1;
9156                                      const double tmp3_1 = w22*y_2;                                      EM_F[INDEX2(k,6,numEq)]+=tmp1 + 6*w5*y_1 + 6*w6*y_2;
9157                                      const double tmp5_1 = w20*y_1;                                      EM_F[INDEX2(k,7,numEq)]+=tmp0 + 6*w5*y_0 + 6*w6*y_3;
9158                                      const double tmp9_1 = w20*y_3;                                  }
9159                                      const double tmp4_1 = tmp1_0*w21;                              } else { // constant data
9160                                      const double tmp8_1 = w22*y_0;                                  for (index_t k=0; k<numEq; k++) {
9161                                      const double tmp2_1 = w20*y_0;                                      EM_F[INDEX2(k,2,numEq)]+=36*w7*y_p[k];
9162                                      const double tmp7_1 = w20*y_2;                                      EM_F[INDEX2(k,3,numEq)]+=36*w7*y_p[k];
9163                                      const double tmp1_1 = tmp0_0*w21;                                      EM_F[INDEX2(k,6,numEq)]+=36*w7*y_p[k];
9164                                      EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_F[INDEX2(k,7,numEq)]+=36*w7*y_p[k];
                                     EM_F[INDEX2(k,3,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_F[INDEX2(k,6,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                     EM_F[INDEX2(k,7,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;  
                                 }  
                             } else { /* constant data */  
                                 for (index_t k=0; k<numEq; k++) {  
                                     const double y_0 = y_p[k];  
                                     const double tmp0_1 = w23*y_0;  
                                     EM_F[INDEX2(k,2,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,6,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,7,numEq)]+=tmp0_1;  
9165                                  }                                  }
9166                              }                              }
9167                          }                          }
9168                            /* GENERATOR SNIP_PDEBC_SYSTEM_3 BOTTOM */
9169                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*(m_NN[1]-2)+k0;                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+m_NN[0]*(m_NN[1]-2)+k0;
9170                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9171                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 9518  void Brick::assemblePDEBoundarySystem(Pa Line 9182  void Brick::assemblePDEBoundarySystem(Pa
9182                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
9183                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
9184                          const index_t e = m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]);                          const index_t e = m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]);
9185                            /* GENERATOR SNIP_PDEBC_SYSTEM_4 TOP */
9186                          ///////////////                          ///////////////
9187                          // process d //                          // process d //
9188                          ///////////////                          ///////////////
# Line 9526  void Brick::assemblePDEBoundarySystem(Pa Line 9191  void Brick::assemblePDEBoundarySystem(Pa
9191                              if (d.actsExpanded()) {                              if (d.actsExpanded()) {
9192                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9193                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9194                                          const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                          const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
9195                                          const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                          const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
9196                                          const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];                                          const double d_2 = d_p[INDEX3(k,m,2,numEq,numComp)];
9197                                          const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];                                          const double d_3 = d_p[INDEX3(k,m,3,numEq,numComp)];
9198                                          const double tmp0_0 = d_1 + d_3;                                          const double tmp0 = w10*(d_0 + d_2);
9199                                          const double tmp2_0 = d_0 + d_1 + d_2 + d_3;                                          const double tmp1 = w11*(d_1 + d_3);
9200                                          const double tmp1_0 = d_0 + d_2;                                          const double tmp2 = w12*(d_0 + d_1 + d_2 + d_3);
9201                                          const double tmp6_0 = d_0 + d_1;                                          const double tmp3 = w12*(d_1 + d_2);
9202                                          const double tmp4_0 = d_0 + d_3;                                          const double tmp4 = w10*(d_1 + d_3);
9203                                          const double tmp5_0 = d_2 + d_3;                                          const double tmp5 = w11*(d_0 + d_2);
9204                                          const double tmp3_0 = d_1 + d_2;                                          const double tmp6 = w12*(d_0 + d_3);
9205                                          const double tmp18_1 = tmp5_0*w24;                                          const double tmp7 = w10*(d_0 + d_1);
9206                                          const double tmp6_1 = tmp1_0*w25;                                          const double tmp8 = w11*(d_2 + d_3);
9207                                          const double tmp4_1 = d_0*w27;                                          const double tmp9 = w10*(d_2 + d_3);
9208                                          const double tmp12_1 = d_2*w27;                                          const double tmp10 = w11*(d_0 + d_1);
9209                                          const double tmp0_1 = tmp0_0*w25;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=d_0*w14 + d_3*w13 + tmp3;
9210                                          const double tmp5_1 = tmp3_0*w26;                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp10 + tmp9;
9211                                          const double tmp2_1 = tmp2_0*w26;                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp4 + tmp5;
9212                                          const double tmp17_1 = tmp6_0*w25;                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp2;
9213                                          const double tmp14_1 = tmp6_0*w24;                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp10 + tmp9;
9214                                          const double tmp11_1 = d_1*w28;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=d_1*w14 + d_2*w13 + tmp6;
9215                                          const double tmp9_1 = d_1*w27;                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp2;
9216                                          const double tmp16_1 = d_3*w27;                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0 + tmp1;
9217                                          const double tmp8_1 = d_2*w28;                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp4 + tmp5;
9218                                          const double tmp7_1 = tmp0_0*w24;                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp2;
9219                                          const double tmp15_1 = d_0*w28;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=d_1*w13 + d_2*w14 + tmp6;
9220                                          const double tmp13_1 = tmp5_0*w25;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp7 + tmp8;
9221                                          const double tmp3_1 = d_3*w28;                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp2;
9222                                          const double tmp10_1 = tmp4_0*w26;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0 + tmp1;
9223                                          const double tmp1_1 = tmp1_0*w24;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp7 + tmp8;
9224                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=d_0*w13 + d_3*w14 + tmp3;
                                         EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp2_1;  
                                         EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp2_1;  
                                         EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp2_1;  
                                         EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                         EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp6_1 + tmp7_1;  
                                         EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp10_1 + tmp8_1 + tmp9_1;  
                                         EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp12_1;  
                                         EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp13_1 + tmp14_1;  
                                         EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp15_1 + tmp16_1 + tmp5_1;  
                                         EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp13_1 + tmp14_1;  
                                         EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp2_1;  
                                         EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp17_1 + tmp18_1;  
                                         EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp17_1 + tmp18_1;  
                                         EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1 + tmp1_1;  
                                         EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp6_1 + tmp7_1;  
9225                                      }                                      }
9226                                   }                                   }
9227                              } else { /* constant data */                              } else { // constant data
9228                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9229                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9230                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double d_0 = d_p[INDEX2(k, m, numEq)];
9231                                          const double tmp2_1 = d_0*w31;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=16*d_0*w12;
9232                                          const double tmp1_1 = d_0*w30;                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=8*d_0*w12;
9233                                          const double tmp0_1 = d_0*w29;                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=8*d_0*w12;
9234                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=4*d_0*w12;
9235                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=8*d_0*w12;
9236                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=16*d_0*w12;
9237                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=4*d_0*w12;
9238                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=8*d_0*w12;
9239                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=8*d_0*w12;
9240                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=4*d_0*w12;
9241                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=16*d_0*w12;
9242                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=8*d_0*w12;
9243                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=4*d_0*w12;
9244                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=8*d_0*w12;
9245                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=8*d_0*w12;
9246                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=16*d_0*w12;
                                         EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1;  
9247                                      }                                      }
9248                                  }                                  }
9249                              }                              }
# Line 9612  void Brick::assemblePDEBoundarySystem(Pa Line 9259  void Brick::assemblePDEBoundarySystem(Pa
9259                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];
9260                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];
9261                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];
9262                                      const double tmp0_0 = y_1 + y_2;                                      const double tmp0 = 6*w12*(y_1 + y_2);
9263                                      const double tmp1_0 = y_0 + y_3;                                      const double tmp1 = 6*w12*(y_0 + y_3);
9264                                      const double tmp7_1 = w34*y_1;                                      EM_F[INDEX2(k,0,numEq)]+=tmp0 + 6*w10*y_3 + 6*w11*y_0;
9265                                      const double tmp3_1 = w32*y_1;                                      EM_F[INDEX2(k,1,numEq)]+=tmp1 + 6*w10*y_2 + 6*w11*y_1;
9266                                      const double tmp8_1 = w32*y_3;                                      EM_F[INDEX2(k,2,numEq)]+=tmp1 + 6*w10*y_1 + 6*w11*y_2;
9267                                      const double tmp0_1 = w32*y_0;                                      EM_F[INDEX2(k,3,numEq)]+=tmp0 + 6*w10*y_0 + 6*w11*y_3;
9268                                      const double tmp2_1 = w34*y_3;                                  }
9269                                      const double tmp9_1 = w34*y_0;                              } else { // constant data
9270                                      const double tmp6_1 = w32*y_2;                                  for (index_t k=0; k<numEq; k++) {
9271                                      const double tmp5_1 = w34*y_2;                                      EM_F[INDEX2(k,0,numEq)]+=36*w12*y_p[k];
9272                                      const double tmp1_1 = tmp0_0*w33;                                      EM_F[INDEX2(k,1,numEq)]+=36*w12*y_p[k];
9273                                      const double tmp4_1 = tmp1_0*w33;                                      EM_F[INDEX2(k,2,numEq)]+=36*w12*y_p[k];
9274                                      EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_F[INDEX2(k,3,numEq)]+=36*w12*y_p[k];
                                     EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                     EM_F[INDEX2(k,3,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;  
                                 }  
                             } else { /* constant data */  
                                 for (index_t k=0; k<numEq; k++) {  
                                     const double y_0 = y_p[k];  
                                     const double tmp0_1 = w35*y_0;  
                                     EM_F[INDEX2(k,0,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,1,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,2,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
9275                                  }                                  }
9276                              }                              }
9277                          }                          }
9278                            /* GENERATOR SNIP_PDEBC_SYSTEM_4 BOTTOM */
9279                          const index_t firstNode=m_NN[0]*k1+k0;                          const index_t firstNode=m_NN[0]*k1+k0;
9280                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9281                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 9656  void Brick::assemblePDEBoundarySystem(Pa Line 9292  void Brick::assemblePDEBoundarySystem(Pa
9292                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
9293                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
9294                          const index_t e = m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]);                          const index_t e = m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]);
9295                            /* GENERATOR SNIP_PDEBC_SYSTEM_5 TOP */
9296                          ///////////////                          ///////////////
9297                          // process d //                          // process d //
9298                          ///////////////                          ///////////////
# Line 9664  void Brick::assemblePDEBoundarySystem(Pa Line 9301  void Brick::assemblePDEBoundarySystem(Pa
9301                              if (d.actsExpanded()) {                              if (d.actsExpanded()) {
9302                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9303                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9304                                          const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                          const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
9305                                          const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                          const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
9306                                          const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];                                          const double d_2 = d_p[INDEX3(k,m,2,numEq,numComp)];
9307                                          const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];                                          const double d_3 = d_p[INDEX3(k,m,3,numEq,numComp)];
9308                                          const double tmp2_0 = d_1 + d_3;                                          const double tmp0 = w12*(d_0 + d_1 + d_2 + d_3);
9309                                          const double tmp0_0 = d_0 + d_1 + d_2 + d_3;                                          const double tmp1 = w10*(d_1 + d_3);
9310                                          const double tmp1_0 = d_0 + d_2;                                          const double tmp2 = w11*(d_0 + d_2);
9311                                          const double tmp3_0 = d_0 + d_1;                                          const double tmp3 = w10*(d_2 + d_3);
9312                                          const double tmp6_0 = d_0 + d_3;                                          const double tmp4 = w11*(d_0 + d_1);
9313                                          const double tmp4_0 = d_2 + d_3;                                          const double tmp5 = w10*(d_0 + d_1);
9314                                          const double tmp5_0 = d_1 + d_2;                                          const double tmp6 = w11*(d_2 + d_3);
9315                                          const double tmp1_1 = tmp1_0*w25;                                          const double tmp7 = w12*(d_1 + d_2);
9316                                          const double tmp11_1 = d_0*w27;                                          const double tmp8 = w10*(d_0 + d_2);
9317                                          const double tmp9_1 = tmp5_0*w26;                                          const double tmp9 = w11*(d_1 + d_3);
9318                                          const double tmp12_1 = tmp2_0*w25;                                          const double tmp10 = w12*(d_0 + d_3);
9319                                          const double tmp3_1 = tmp3_0*w25;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=d_0*w14 + d_3*w13 + tmp7;
9320                                          const double tmp6_1 = tmp3_0*w24;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp3 + tmp4;
9321                                          const double tmp2_1 = tmp2_0*w24;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp1 + tmp2;
9322                                          const double tmp10_1 = d_3*w28;                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0;
9323                                          const double tmp16_1 = tmp6_0*w26;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp3 + tmp4;
9324                                          const double tmp17_1 = d_1*w28;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=d_1*w14 + d_2*w13 + tmp10;
9325                                          const double tmp7_1 = d_0*w28;                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0;
9326                                          const double tmp14_1 = d_2*w28;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp8 + tmp9;
9327                                          const double tmp18_1 = d_2*w27;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp1 + tmp2;
9328                                          const double tmp15_1 = d_1*w27;                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0;
9329                                          const double tmp0_1 = tmp0_0*w26;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=d_1*w13 + d_2*w14 + tmp10;
9330                                          const double tmp4_1 = tmp4_0*w24;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp5 + tmp6;
9331                                          const double tmp5_1 = tmp4_0*w25;                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0;
9332                                          const double tmp8_1 = d_3*w27;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp8 + tmp9;
9333                                          const double tmp13_1 = tmp1_0*w24;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp5 + tmp6;
9334                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=d_0*w13 + d_3*w14 + tmp7;
                                         EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp1_1 + tmp2_1;  
                                         EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1;  
                                         EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp5_1 + tmp6_1;  
                                         EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp5_1 + tmp6_1;  
                                         EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                         EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp9_1;  
                                         EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp12_1 + tmp13_1;  
                                         EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp14_1 + tmp15_1 + tmp16_1;  
                                         EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp3_1 + tmp4_1;  
                                         EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp16_1 + tmp17_1 + tmp18_1;  
                                         EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp12_1 + tmp13_1;  
                                         EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp1_1 + tmp2_1;  
                                         EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0_1;  
9335                                      }                                      }
9336                                   }                                   }
9337                              } else { /* constant data */                              } else { // constant data
9338                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9339                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9340                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double d_0 = d_p[INDEX2(k, m, numEq)];
9341                                          const double tmp2_1 = d_0*w31;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=16*d_0*w12;
9342                                          const double tmp0_1 = d_0*w30;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=8*d_0*w12;
9343                                          const double tmp1_1 = d_0*w29;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=8*d_0*w12;
9344                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=4*d_0*w12;
9345                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=8*d_0*w12;
9346                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=16*d_0*w12;
9347                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1;                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=4*d_0*w12;
9348                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=8*d_0*w12;
9349                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=8*d_0*w12;
9350                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=4*d_0*w12;
9351                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=16*d_0*w12;
9352                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=8*d_0*w12;
9353                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=4*d_0*w12;
9354                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=8*d_0*w12;
9355                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp2_1;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=8*d_0*w12;
9356                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp1_1;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=16*d_0*w12;
                                         EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0_1;  
                                         EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp1_1;  
                                         EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0_1;  
9357                                      }                                      }
9358                                  }                                  }
9359                              }                              }
# Line 9750  void Brick::assemblePDEBoundarySystem(Pa Line 9369  void Brick::assemblePDEBoundarySystem(Pa
9369                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];                                      const double y_1 = y_p[INDEX2(k, 1, numEq)];
9370                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];                                      const double y_2 = y_p[INDEX2(k, 2, numEq)];
9371                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];                                      const double y_3 = y_p[INDEX2(k, 3, numEq)];
9372                                      const double tmp0_0 = y_1 + y_2;                                      const double tmp0 = 6*w12*(y_1 + y_2);
9373                                      const double tmp1_0 = y_0 + y_3;                                      const double tmp1 = 6*w12*(y_0 + y_3);
9374                                      const double tmp7_1 = w34*y_1;                                      EM_F[INDEX2(k,4,numEq)]+=tmp0 + 6*w10*y_3 + 6*w11*y_0;
9375                                      const double tmp3_1 = w32*y_1;                                      EM_F[INDEX2(k,5,numEq)]+=tmp1 + 6*w10*y_2 + 6*w11*y_1;
9376                                      const double tmp8_1 = w32*y_3;                                      EM_F[INDEX2(k,6,numEq)]+=tmp1 + 6*w10*y_1 + 6*w11*y_2;
9377                                      const double tmp0_1 = w32*y_0;                                      EM_F[INDEX2(k,7,numEq)]+=tmp0 + 6*w10*y_0 + 6*w11*y_3;
9378                                      const double tmp2_1 = w34*y_3;                                  }
9379                                      const double tmp9_1 = w34*y_0;                              } else { // constant data
9380                                      const double tmp6_1 = w32*y_2;                                  for (index_t k=0; k<numEq; k++) {
9381                                      const double tmp5_1 = w34*y_2;                                      EM_F[INDEX2(k,4,numEq)]+=36*w12*y_p[k];
9382                                      const double tmp1_1 = tmp0_0*w33;                                      EM_F[INDEX2(k,5,numEq)]+=36*w12*y_p[k];
9383                                      const double tmp4_1 = tmp1_0*w33;                                      EM_F[INDEX2(k,6,numEq)]+=36*w12*y_p[k];
9384                                      EM_F[INDEX2(k,4,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                      EM_F[INDEX2(k,7,numEq)]+=36*w12*y_p[k];
                                     EM_F[INDEX2(k,5,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_F[INDEX2(k,6,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                     EM_F[INDEX2(k,7,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;  
                                 }  
                             } else { /* constant data */  
                                 for (index_t k=0; k<numEq; k++) {  
                                     const double y_0 = y_p[k];  
                                     const double tmp0_1 = w35*y_0;  
                                     EM_F[INDEX2(k,4,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,5,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,6,numEq)]+=tmp0_1;  
                                     EM_F[INDEX2(k,7,numEq)]+=tmp0_1;  
9385                                  }                                  }
9386                              }                              }
9387                          }                          }
9388                            /* GENERATOR SNIP_PDEBC_SYSTEM_5 BOTTOM */
9389                          const index_t firstNode=m_NN[0]*m_NN[1]*(m_NN[2]-2)+m_NN[0]*k1+k0;                          const index_t firstNode=m_NN[0]*m_NN[1]*(m_NN[2]-2)+m_NN[0]*k1+k0;
9390                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9391                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 9799  void Brick::assemblePDEBoundarySystemRed Line 9407  void Brick::assemblePDEBoundarySystemRed
9407          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
9408          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
9409      }      }
9410      const double w0 = 0.0625*m_dx[0]*m_dx[1];      const double w0 = m_dx[0]*m_dx[1]/16.;
9411      const double w1 = 0.0625*m_dx[0]*m_dx[2];      const double w1 = m_dx[0]*m_dx[2]/16.;
9412      const double w2 = 0.0625*m_dx[1]*m_dx[2];      const double w2 = m_dx[1]*m_dx[2]/16.;
9413      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
9414      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
9415      rhs.requireWrite();      rhs.requireWrite();
# Line 9823  void Brick::assemblePDEBoundarySystemRed Line 9431  void Brick::assemblePDEBoundarySystemRed
9431                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
9432                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
9433                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w2;                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w2;
                                     EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0;  
9434                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0;
9435                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0;
9436                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0;
9437                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp0;
                                     EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0;  
9438                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0;
9439                                      EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0;
9440                                        EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp0;
9441                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0;
9442                                        EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0;
9443                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp0;
9444                                        EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0;
9445                                        EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0;
9446                                        EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp0;
9447                                        EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0;
9448                                        EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0;
9449                                        EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0;
9450                                  }                                  }
9451                              }                              }
9452                          }                          }
# Line 9878  void Brick::assemblePDEBoundarySystemRed Line 9486  void Brick::assemblePDEBoundarySystemRed
9486                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
9487                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
9488                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w2;                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w2;
9489                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0;
9490                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0;
9491                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0;
9492                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp0;
9493                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0;
9494                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0;
                                     EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0;  
9495                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp0;
9496                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0;
9497                                      EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0;
9498                                        EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp0;
9499                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0;
9500                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0;
                                     EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0;  
9501                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp0;
9502                                        EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0;
9503                                        EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0;
9504                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0;
9505                                  }                                  }
9506                              }                              }
9507                          }                          }
# Line 9933  void Brick::assemblePDEBoundarySystemRed Line 9541  void Brick::assemblePDEBoundarySystemRed
9541                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
9542                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
9543                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w1;                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w1;
                                     EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0;  
9544                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0;
9545                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0;
9546                                        EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0;
9547                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp0;
9548                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0;
9549                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0;
9550                                        EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp0;
9551                                        EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0;
9552                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0;
                                     EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0;  
9553                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp0;
9554                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0;
9555                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0;
9556                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp0;
9557                                        EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0;
9558                                        EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0;
9559                                        EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0;
9560                                  }                                  }
9561                              }                              }
9562                          }                          }
# Line 9988  void Brick::assemblePDEBoundarySystemRed Line 9596  void Brick::assemblePDEBoundarySystemRed
9596                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
9597                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
9598                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w1;                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w1;
9599                                      EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0;
9600                                      EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0;
9601                                      EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0;
9602                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0;
9603                                      EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0;
9604                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0;
                                     EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0;  
9605                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0;
9606                                        EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0;
9607                                        EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0;
9608                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0;
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0;  
9609                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0;
9610                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0;
9611                                      EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0;
9612                                        EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0;
9613                                        EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0;
9614                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0;
9615                                  }                                  }
9616                              }                              }
9617                          }                          }
# Line 10043  void Brick::assemblePDEBoundarySystemRed Line 9651  void Brick::assemblePDEBoundarySystemRed
9651                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
9652                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
9653                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w0;                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w0;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0;  
9654                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0;
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0;  
9655                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0;
9656                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0;
9657                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp0;
9658                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0;
9659                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0;
9660                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0;
9661                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0;
9662                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0;
9663                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0;
9664                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0;
9665                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0;
9666                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp0;
9667                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0;
9668                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0;
9669                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0;
9670                                  }                                  }
9671                              }                              }
9672                          }                          }
# Line 10098  void Brick::assemblePDEBoundarySystemRed Line 9706  void Brick::assemblePDEBoundarySystemRed
9706                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
9707                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
9708                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w0;                                      const double tmp0 = d_p[INDEX2(k, m, numEq)]*w0;
                                     EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0;  
                                     EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0;  
9709                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0;
9710                                      EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0;
9711                                      EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0;
9712                                      EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0;
9713                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0;
9714                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0;
                                     EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0;  
9715                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0;
9716                                        EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0;
9717                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0;
9718                                      EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0;                                      EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0;
9719                                        EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0;
9720                                        EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0;
9721                                        EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0;
9722                                        EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0;
9723                                        EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0;
9724                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0;
9725                                  }                                  }
9726                              }                              }
9727                          }                          }

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

  ViewVC Help
Powered by ViewVC 1.1.26