/[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 4378 by caltinay, Wed Apr 24 06:21:42 2013 UTC revision 4382 by caltinay, Thu Apr 25 04:44:14 2013 UTC
# Line 5264  void Brick::assemblePDESystem(Paso_Syste Line 5264  void Brick::assemblePDESystem(Paso_Syste
5264      const double w1 = w8*(-SQRT3 + 2);      const double w1 = w8*(-SQRT3 + 2);
5265      const double w20 = w8*(4*SQRT3 - 7);      const double w20 = w8*(4*SQRT3 - 7);
5266      const double w21 = w8*(-4*SQRT3 - 7);      const double w21 = w8*(-4*SQRT3 - 7);
   
5267      const double w54 = -m_dx[0]*m_dx[1]/72;      const double w54 = -m_dx[0]*m_dx[1]/72;
5268      const double w68 = -m_dx[0]*m_dx[1]/48;      const double w68 = -m_dx[0]*m_dx[1]/48;
5269      const double w38 = w68*(-SQRT3 - 3)/36;      const double w38 = w68*(-SQRT3 - 3)/36;
# Line 5275  void Brick::assemblePDESystem(Paso_Syste Line 5274  void Brick::assemblePDESystem(Paso_Syste
5274      const double w44 = w68*(19*SQRT3 - 33)/36;      const double w44 = w68*(19*SQRT3 - 33)/36;
5275      const double w66 = w68*(SQRT3 + 2);      const double w66 = w68*(SQRT3 + 2);
5276      const double w70 = w68*(-SQRT3 + 2);      const double w70 = w68*(-SQRT3 + 2);
   
5277      const double w56 = -m_dx[0]*m_dx[2]/72;      const double w56 = -m_dx[0]*m_dx[2]/72;
5278      const double w67 = -m_dx[0]*m_dx[2]/48;      const double w67 = -m_dx[0]*m_dx[2]/48;
5279      const double w37 = w67*(-SQRT3 - 3)/36;      const double w37 = w67*(-SQRT3 - 3)/36;
# Line 5286  void Brick::assemblePDESystem(Paso_Syste Line 5284  void Brick::assemblePDESystem(Paso_Syste
5284      const double w48 = w67*(-19*SQRT3 + 33)/36;      const double w48 = w67*(-19*SQRT3 + 33)/36;
5285      const double w65 = w67*(SQRT3 + 2);      const double w65 = w67*(SQRT3 + 2);
5286      const double w71 = w67*(-SQRT3 + 2);      const double w71 = w67*(-SQRT3 + 2);
   
5287      const double w55 = -m_dx[1]*m_dx[2]/72;      const double w55 = -m_dx[1]*m_dx[2]/72;
5288      const double w69 = -m_dx[1]*m_dx[2]/48;      const double w69 = -m_dx[1]*m_dx[2]/48;
5289      const double w36 = w69*(SQRT3 - 3)/36;      const double w36 = w69*(SQRT3 - 3)/36;
# Line 5297  void Brick::assemblePDESystem(Paso_Syste Line 5294  void Brick::assemblePDESystem(Paso_Syste
5294      const double w50 = w69*(-19*SQRT3 - 33)/36;      const double w50 = w69*(-19*SQRT3 - 33)/36;
5295      const double w64 = w69*(SQRT3 + 2);      const double w64 = w69*(SQRT3 + 2);
5296      const double w72 = w69*(-SQRT3 + 2);      const double w72 = w69*(-SQRT3 + 2);
   
5297      const double w58 = m_dx[0]*m_dx[1]*m_dx[2]/1728;      const double w58 = m_dx[0]*m_dx[1]*m_dx[2]/1728;
5298        const double w60 = w58*(-SQRT3 + 2);
5299      const double w61 = w58*(SQRT3 + 2);      const double w61 = w58*(SQRT3 + 2);
     const double w62 = w58*(15*SQRT3 + 26);  
5300      const double w57 = w58*(-4*SQRT3 + 7);      const double w57 = w58*(-4*SQRT3 + 7);
     const double w60 = w58*(-SQRT3 + 2);  
     const double w63 = w58*(-15*SQRT3 + 26);  
5301      const double w59 = w58*(4*SQRT3 + 7);      const double w59 = w58*(4*SQRT3 + 7);
5302        const double w62 = w58*(15*SQRT3 + 26);
5303        const double w63 = w58*(-15*SQRT3 + 26);
5304      const double w75 = w58*6*(SQRT3 + 3);      const double w75 = w58*6*(SQRT3 + 3);
5305      const double w76 = w58*6*(-SQRT3 + 3);      const double w76 = w58*6*(-SQRT3 + 3);
5306      const double w74 = w58*6*(5*SQRT3 + 9);      const double w74 = w58*6*(5*SQRT3 + 9);
5307      const double w77 = w58*6*(-5*SQRT3 + 9);      const double w77 = w58*6*(-5*SQRT3 + 9);
   
5308      const double w13 = -m_dx[0]*m_dx[1]/(288*m_dx[2]);      const double w13 = -m_dx[0]*m_dx[1]/(288*m_dx[2]);
5309      const double w19 = w13*(4*SQRT3 + 7);      const double w19 = w13*(4*SQRT3 + 7);
5310      const double w7 = w13*(-4*SQRT3 + 7);      const double w7 = w13*(-4*SQRT3 + 7);
     const double w25 = w13*(-SQRT3 - 2);  
5311      const double w23 = w13*(+SQRT3 - 2);      const double w23 = w13*(+SQRT3 - 2);
5312      const double w32 = -8*w13;      const double w25 = w13*(-SQRT3 - 2);
   
5313      const double w22 = -m_dx[0]*m_dx[2]/(288*m_dx[1]);      const double w22 = -m_dx[0]*m_dx[2]/(288*m_dx[1]);
5314      const double w3 = w22*(SQRT3 - 2);      const double w3 = w22*(SQRT3 - 2);
5315      const double w9 = w22*(-SQRT3 - 2);      const double w9 = w22*(-SQRT3 - 2);
5316      const double w24 = w22*(4*SQRT3 + 7);      const double w24 = w22*(4*SQRT3 + 7);
5317      const double w26 = w22*(-4*SQRT3 + 7);      const double w26 = w22*(-4*SQRT3 + 7);
     const double w31 = -8*w22;  
   
5318      const double w27 = -m_dx[1]*m_dx[2]/(288*m_dx[0]);      const double w27 = -m_dx[1]*m_dx[2]/(288*m_dx[0]);
     const double w14 = w27*(-SQRT3 - 2);  
5319      const double w0 = w27*(SQRT3 - 2);      const double w0 = w27*(SQRT3 - 2);
5320        const double w14 = w27*(-SQRT3 - 2);
5321      const double w28 = w27*(-4*SQRT3 + 7);      const double w28 = w27*(-4*SQRT3 + 7);
5322      const double w29 = w27*(4*SQRT3 + 7);      const double w29 = w27*(4*SQRT3 + 7);
     const double w30 = -8*w27;  
5323    
5324      rhs.requireWrite();      rhs.requireWrite();
5325  #pragma omp parallel  #pragma omp parallel
# Line 6004  void Brick::assemblePDESystem(Paso_Syste Line 5994  void Brick::assemblePDESystem(Paso_Syste
5994                              } else { // constant data                              } else { // constant data
5995                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
5996                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
5997                                          const double A_00 = A_p[INDEX4(k,0,m,0, numEq,3, numComp)];                                          const double Aw00 = 8*A_p[INDEX4(k,0,m,0, numEq,3, numComp)]*w27;
5998                                          const double A_01 = A_p[INDEX4(k,0,m,1, numEq,3, numComp)];                                          const double Aw01 = 12*A_p[INDEX4(k,0,m,1, numEq,3, numComp)]*w8;
5999                                          const double A_02 = A_p[INDEX4(k,0,m,2, numEq,3, numComp)];                                          const double Aw02 = 12*A_p[INDEX4(k,0,m,2, numEq,3, numComp)]*w11;
6000                                          const double A_10 = A_p[INDEX4(k,1,m,0, numEq,3, numComp)];                                          const double Aw10 = 12*A_p[INDEX4(k,1,m,0, numEq,3, numComp)]*w8;
6001                                          const double A_11 = A_p[INDEX4(k,1,m,1, numEq,3, numComp)];                                          const double Aw11 = 8*A_p[INDEX4(k,1,m,1, numEq,3, numComp)]*w22;
6002                                          const double A_12 = A_p[INDEX4(k,1,m,2, numEq,3, numComp)];                                          const double Aw12 = 12*A_p[INDEX4(k,1,m,2, numEq,3, numComp)]*w10;
6003                                          const double A_20 = A_p[INDEX4(k,2,m,0, numEq,3, numComp)];                                          const double Aw20 = 12*A_p[INDEX4(k,2,m,0, numEq,3, numComp)]*w11;
6004                                          const double A_21 = A_p[INDEX4(k,2,m,1, numEq,3, numComp)];                                          const double Aw21 = 12*A_p[INDEX4(k,2,m,1, numEq,3, numComp)]*w10;
6005                                          const double A_22 = A_p[INDEX4(k,2,m,2, numEq,3, numComp)];                                          const double Aw22 = 8*A_p[INDEX4(k,2,m,2, numEq,3, numComp)]*w13;
6006                                          const double tmp0 = 24*w11*(-A_02 + A_20);                                          const double tmp0 = Aw01 + Aw10;
6007                                          const double tmp1 = 24*w10*(A_12 - A_21);                                          const double tmp1 = Aw01 - Aw10;
6008                                          const double tmp2 = 12*w8*(A_01 + A_10);                                          const double tmp2 = Aw02 + Aw20;
6009                                          const double tmp3 = 12*w11*(-A_02 + A_20);                                          const double tmp3 = Aw02 - Aw20;
6010                                          const double tmp4 = 12*w10*(A_12 - A_21);                                          const double tmp4 = Aw12 + Aw21;
6011                                          const double tmp5 = 24*w8*(-A_01 - A_10);                                          const double tmp5 = Aw12 - Aw21;
6012                                          const double tmp6 = 12*w11*(-A_02 - A_20);                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=-4*Aw00 - 4*Aw11 - 4*Aw22 + 2*tmp0 + 2*tmp2 - 2*tmp4;
6013                                          const double tmp7 = 24*w10*(-A_12 + A_21);                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+= 4*Aw00 - 2*Aw11 - 2*Aw22 - tmp4 + 2*tmp1 + 2*tmp3;
6014                                          const double tmp8 = 24*w8*(A_01 - A_10);                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=-2*Aw00 + 4*Aw11 - 2*Aw22 - 2*tmp1 + tmp2 - 2*tmp5;
6015                                          const double tmp9 = 24*w11*(A_02 - A_20);                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+= 2*Aw00 + 2*Aw11 -   Aw22 - 2*tmp0 + tmp3 - tmp5;
6016                                          const double tmp10 = 12*w10*(A_12 + A_21);                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=-2*Aw00 - 2*Aw11 + 4*Aw22 - 2*tmp3 + 2*tmp5 + tmp0;
6017                                          const double tmp11 = 24*w8*(-A_01 + A_10);                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+= 2*Aw00 -   Aw11 + 2*Aw22 - 2*tmp2 + tmp1 + tmp5;
6018                                          const double tmp12 = 12*w8*(-A_01 - A_10);                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=  -Aw00 + 2*Aw11 + 2*Aw22 + 2*tmp4 - tmp1 - tmp3;
6019                                          const double tmp13 = 12*w11*(A_02 - A_20);                                          EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=   Aw00 +   Aw11 +   Aw22 + tmp4 - tmp0 - tmp2;
6020                                          const double tmp14 = 24*w8*(A_01 + A_10);                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+= 4*Aw00 - 2*Aw11 - 2*Aw22 - 2*tmp3 - 2*tmp1 - tmp4;
6021                                          const double tmp15 = 12*w11*(A_02 + A_20);                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=-4*Aw00 - 4*Aw11 - 4*Aw22 - 2*tmp2 - 2*tmp4 - 2*tmp0;
6022                                          const double tmp16 = 12*w10*(-A_12 - A_21);                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+= 2*Aw00 + 2*Aw11 -   Aw22 + 2*tmp0 - tmp5 - tmp3;
6023                                          const double tmp17 = 12*w10*(-A_12 + A_21);                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=-2*Aw00 + 4*Aw11 - 2*Aw22 - tmp2 - 2*tmp5 + 2*tmp1;
6024                                          const double tmp18 = 24*w11*(-A_02 - A_20);                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+= 2*Aw00 -   Aw11 + 2*Aw22 + 2*tmp2 - tmp1 + tmp5;
6025                                          const double tmp19 = 12*w8*(A_01 - A_10);                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=-2*Aw00 - 2*Aw11 + 4*Aw22 + 2*tmp5 - tmp0 + 2*tmp3;
6026                                          const double tmp20 = 24*w10*(A_12 + A_21);                                          EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=   Aw00 +   Aw11 +   Aw22 + tmp4 + tmp2 + tmp0;
6027                                          const double tmp21 = 24*w11*(A_02 + A_20);                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=  -Aw00 + 2*Aw11 + 2*Aw22 + tmp3 + tmp1 + 2*tmp4;
6028                                          const double tmp22 = 12*w8*(-A_01 + A_10);                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=-2*Aw00 + 4*Aw11 - 2*Aw22 + 2*tmp5 + tmp2 + 2*tmp1;
6029                                          const double tmp23 = 24*w10*(-A_12 - A_21);                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+= 2*Aw00 + 2*Aw11 -   Aw22 + tmp3 + 2*tmp0 + tmp5;
6030                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=4*A_00*w30 + 4*A_11*w31 + 4*A_22*w32 + tmp14 + tmp21 + tmp23;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=-4*Aw00 - 4*Aw11 - 4*Aw22 + 2*tmp4 + 2*tmp2 - 2*tmp0;
6031                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=-4*A_00*w30 + 2*A_11*w31 + 2*A_22*w32 + tmp16 + tmp8 + tmp9;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+= 4*Aw00 - 2*Aw11 - 2*Aw22 + tmp4 - 2*tmp1 + 2*tmp3;
6032                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=2*A_00*w30 + 32*A_11*w22 + 2*A_22*w32 + tmp11 + tmp15 + tmp7;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=  -Aw00 + 2*Aw11 + 2*Aw22 + tmp1 - 2*tmp4 - tmp3;
6033                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=-2*A_00*w30 - 2*A_11*w31 + A_22*w32 + tmp13 + tmp17 + tmp5;                                          EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=   Aw00 +   Aw11 +   Aw22 - tmp4 + tmp0 - tmp2;
6034                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=2*A_00*w30 + 2*A_11*w31 + 32*A_22*w13 + tmp0 + tmp1 + tmp2;                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=-2*Aw00 - 2*Aw11 + 4*Aw22 - 2*tmp3 - tmp0 - 2*tmp5;
6035                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=-2*A_00*w30 + A_11*w31 + 16*A_22*w13 + tmp18 + tmp19 + tmp4;                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+= 2*Aw00 -   Aw11 + 2*Aw22 - tmp5 - 2*tmp2 - tmp1;
6036                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=A_00*w30 - 2*A_11*w31 + 16*A_22*w13 + tmp20 + tmp22 + tmp3;                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+= 2*Aw00 + 2*Aw11 -   Aw22 - tmp3 + tmp5 - 2*tmp0;
6037                                          EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=8*A_00*w27 + 8*A_11*w22 - A_22*w32 + tmp10 + tmp12 + tmp6;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=-2*Aw00 + 4*Aw11 - 2*Aw22 + 2*tmp5 - 2*tmp1 - tmp2;
6038                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=-4*A_00*w30 + 2*A_11*w31 + 2*A_22*w32 + tmp0 + tmp11 + tmp16;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+= 4*Aw00 - 2*Aw11 - 2*Aw22 - 2*tmp3 + tmp4 + 2*tmp1;
6039                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=4*A_00*w30 + 4*A_11*w31 + 4*A_22*w32 + tmp18 + tmp23 + tmp5;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=-4*Aw00 - 4*Aw11 - 4*Aw22 + 2*tmp0 - 2*tmp2 + 2*tmp4;
6040                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=-2*A_00*w30 - 2*A_11*w31 + A_22*w32 + tmp14 + tmp17 + tmp3;                                          EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=   Aw00 +   Aw11 +   Aw22 - tmp0 + tmp2 - tmp4;
6041                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=2*A_00*w30 + 32*A_11*w22 + 2*A_22*w32 + tmp6 + tmp7 + tmp8;                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=  -Aw00 + 2*Aw11 + 2*Aw22 + tmp3 - tmp1 - 2*tmp4;
6042                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=-2*A_00*w30 + A_11*w31 + 16*A_22*w13 + tmp21 + tmp22 + tmp4;                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+= 2*Aw00 -   Aw11 + 2*Aw22 - tmp5 + tmp1 + 2*tmp2;
6043                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=2*A_00*w30 + 2*A_11*w31 + 32*A_22*w13 + tmp1 + tmp12 + tmp9;                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=-2*Aw00 - 2*Aw11 + 4*Aw22 + tmp0 - 2*tmp5 + 2*tmp3;
6044                                          EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=8*A_00*w27 + 8*A_11*w22 - A_22*w32 + tmp10 + tmp15 + tmp2;                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=-2*Aw00 - 2*Aw11 + 4*Aw22 + tmp0 - 2*tmp5 + 2*tmp3;
6045                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=A_00*w30 - 2*A_11*w31 + 16*A_22*w13 + tmp13 + tmp19 + tmp20;                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+= 2*Aw00 -   Aw11 + 2*Aw22 - tmp5 + tmp1 + 2*tmp2;
6046                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=2*A_00*w30 + 32*A_11*w22 + 2*A_22*w32 + tmp1 + tmp15 + tmp8;                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=  -Aw00 + 2*Aw11 + 2*Aw22 + tmp3 - tmp1 - 2*tmp4;
6047                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=-2*A_00*w30 - 2*A_11*w31 + A_22*w32 + tmp13 + tmp14 + tmp4;                                          EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=   Aw00 +   Aw11 +   Aw22 - tmp0 + tmp2 - tmp4;
6048                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=4*A_00*w30 + 4*A_11*w31 + 4*A_22*w32 + tmp20 + tmp21 + tmp5;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=-4*Aw00 - 4*Aw11 - 4*Aw22 + 2*tmp0 - 2*tmp2 + 2*tmp4;
6049                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=-4*A_00*w30 + 2*A_11*w31 + 2*A_22*w32 + tmp10 + tmp11 + tmp9;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+= 4*Aw00 - 2*Aw11 - 2*Aw22 - 2*tmp3 + tmp4 + 2*tmp1;
6050                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=A_00*w30 - 2*A_11*w31 + 16*A_22*w13 + tmp19 + tmp23 + tmp3;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=-2*Aw00 + 4*Aw11 - 2*Aw22 + 2*tmp5 - 2*tmp1 - tmp2;
6051                                          EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=8*A_00*w27 + 8*A_11*w22 - A_22*w32 + tmp16 + tmp2 + tmp6;                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+= 2*Aw00 + 2*Aw11 -   Aw22 - tmp3 + tmp5 - 2*tmp0;
6052                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=2*A_00*w30 + 2*A_11*w31 + 32*A_22*w13 + tmp0 + tmp12 + tmp7;                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+= 2*Aw00 -   Aw11 + 2*Aw22 - tmp5 - 2*tmp2 - tmp1;
6053                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=-2*A_00*w30 + A_11*w31 + 16*A_22*w13 + tmp17 + tmp18 + tmp22;                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=-2*Aw00 - 2*Aw11 + 4*Aw22 - 2*tmp3 - tmp0 - 2*tmp5;
6054                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=-2*A_00*w30 - 2*A_11*w31 + A_22*w32 + tmp3 + tmp4 + tmp5;                                          EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=   Aw00 +   Aw11 +   Aw22 - tmp4 + tmp0 - tmp2;
6055                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=2*A_00*w30 + 32*A_11*w22 + 2*A_22*w32 + tmp1 + tmp11 + tmp6;                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=  -Aw00 + 2*Aw11 + 2*Aw22 + tmp1 - 2*tmp4 - tmp3;
6056                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=-4*A_00*w30 + 2*A_11*w31 + 2*A_22*w32 + tmp0 + tmp10 + tmp8;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+= 4*Aw00 - 2*Aw11 - 2*Aw22 + tmp4 - 2*tmp1 + 2*tmp3;
6057                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=4*A_00*w30 + 4*A_11*w31 + 4*A_22*w32 + tmp14 + tmp18 + tmp20;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=-4*Aw00 - 4*Aw11 - 4*Aw22 + 2*tmp4 + 2*tmp2 - 2*tmp0;
6058                                          EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=8*A_00*w27 + 8*A_11*w22 - A_22*w32 + tmp12 + tmp15 + tmp16;                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+= 2*Aw00 + 2*Aw11 -   Aw22 + tmp3 + 2*tmp0 + tmp5;
6059                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=A_00*w30 - 2*A_11*w31 + 16*A_22*w13 + tmp13 + tmp22 + tmp23;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=-2*Aw00 + 4*Aw11 - 2*Aw22 + 2*tmp5 + tmp2 + 2*tmp1;
6060                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=-2*A_00*w30 + A_11*w31 + 16*A_22*w13 + tmp17 + tmp19 + tmp21;                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=  -Aw00 + 2*Aw11 + 2*Aw22 + tmp3 + tmp1 + 2*tmp4;
6061                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=2*A_00*w30 + 2*A_11*w31 + 32*A_22*w13 + tmp2 + tmp7 + tmp9;                                          EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=   Aw00 +   Aw11 +   Aw22 + tmp4 + tmp2 + tmp0;
6062                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=2*A_00*w30 + 2*A_11*w31 + 32*A_22*w13 + tmp2 + tmp7 + tmp9;                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=-2*Aw00 - 2*Aw11 + 4*Aw22 + 2*tmp5 - tmp0 + 2*tmp3;
6063                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=-2*A_00*w30 + A_11*w31 + 16*A_22*w13 + tmp17 + tmp19 + tmp21;                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+= 2*Aw00 -   Aw11 + 2*Aw22 + 2*tmp2 - tmp1 + tmp5;
6064                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=A_00*w30 - 2*A_11*w31 + 16*A_22*w13 + tmp13 + tmp22 + tmp23;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=-2*Aw00 + 4*Aw11 - 2*Aw22 - tmp2 - 2*tmp5 + 2*tmp1;
6065                                          EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=8*A_00*w27 + 8*A_11*w22 - A_22*w32 + tmp12 + tmp15 + tmp16;                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+= 2*Aw00 + 2*Aw11 -   Aw22 + 2*tmp0 - tmp5 - tmp3;
6066                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=4*A_00*w30 + 4*A_11*w31 + 4*A_22*w32 + tmp14 + tmp18 + tmp20;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=-4*Aw00 - 4*Aw11 - 4*Aw22 - 2*tmp2 - 2*tmp4 - 2*tmp0;
6067                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=-4*A_00*w30 + 2*A_11*w31 + 2*A_22*w32 + tmp0 + tmp10 + tmp8;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+= 4*Aw00 - 2*Aw11 - 2*Aw22 - 2*tmp3 - 2*tmp1 - tmp4;
6068                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=2*A_00*w30 + 32*A_11*w22 + 2*A_22*w32 + tmp1 + tmp11 + tmp6;                                          EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=   Aw00 +   Aw11 +   Aw22 + tmp4 - tmp0 - tmp2;
6069                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=-2*A_00*w30 - 2*A_11*w31 + A_22*w32 + tmp3 + tmp4 + tmp5;                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=  -Aw00 + 2*Aw11 + 2*Aw22 + 2*tmp4 - tmp1 - tmp3;
6070                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=-2*A_00*w30 + A_11*w31 + 16*A_22*w13 + tmp17 + tmp18 + tmp22;                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+= 2*Aw00 -   Aw11 + 2*Aw22 - 2*tmp2 + tmp1 + tmp5;
6071                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=2*A_00*w30 + 2*A_11*w31 + 32*A_22*w13 + tmp0 + tmp12 + tmp7;                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=-2*Aw00 - 2*Aw11 + 4*Aw22 - 2*tmp3 + 2*tmp5 + tmp0;
6072                                          EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=8*A_00*w27 + 8*A_11*w22 - A_22*w32 + tmp16 + tmp2 + tmp6;                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+= 2*Aw00 + 2*Aw11 -   Aw22 + tmp3 - tmp5 - 2*tmp0;
6073                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=A_00*w30 - 2*A_11*w31 + 16*A_22*w13 + tmp19 + tmp23 + tmp3;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=-2*Aw00 + 4*Aw11 - 2*Aw22 - 2*tmp1 + tmp2 - 2*tmp5;
6074                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=-4*A_00*w30 + 2*A_11*w31 + 2*A_22*w32 + tmp10 + tmp11 + tmp9;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+= 4*Aw00 - 2*Aw11 - 2*Aw22 - tmp4 + 2*tmp1 + 2*tmp3;
6075                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=4*A_00*w30 + 4*A_11*w31 + 4*A_22*w32 + tmp20 + tmp21 + tmp5;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=-4*Aw00 - 4*Aw11 - 4*Aw22 + 2*tmp0 + 2*tmp2 - 2*tmp4;
                                         EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=-2*A_00*w30 - 2*A_11*w31 + A_22*w32 + tmp13 + tmp14 + tmp4;  
                                         EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=2*A_00*w30 + 32*A_11*w22 + 2*A_22*w32 + tmp1 + tmp15 + tmp8;  
                                         EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=A_00*w30 - 2*A_11*w31 + 16*A_22*w13 + tmp13 + tmp19 + tmp20;  
                                         EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=8*A_00*w27 + 8*A_11*w22 - A_22*w32 + tmp10 + tmp15 + tmp2;  
                                         EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=2*A_00*w30 + 2*A_11*w31 + 32*A_22*w13 + tmp1 + tmp12 + tmp9;  
                                         EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=-2*A_00*w30 + A_11*w31 + 16*A_22*w13 + tmp21 + tmp22 + tmp4;  
                                         EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=2*A_00*w30 + 32*A_11*w22 + 2*A_22*w32 + tmp6 + tmp7 + tmp8;  
                                         EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=-2*A_00*w30 - 2*A_11*w31 + A_22*w32 + tmp14 + tmp17 + tmp3;  
                                         EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=4*A_00*w30 + 4*A_11*w31 + 4*A_22*w32 + tmp18 + tmp23 + tmp5;  
                                         EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=-4*A_00*w30 + 2*A_11*w31 + 2*A_22*w32 + tmp0 + tmp11 + tmp16;  
                                         EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=8*A_00*w27 + 8*A_11*w22 - A_22*w32 + tmp10 + tmp12 + tmp6;  
                                         EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=A_00*w30 - 2*A_11*w31 + 16*A_22*w13 + tmp20 + tmp22 + tmp3;  
                                         EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=-2*A_00*w30 + A_11*w31 + 16*A_22*w13 + tmp18 + tmp19 + tmp4;  
                                         EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=2*A_00*w30 + 2*A_11*w31 + 32*A_22*w13 + tmp0 + tmp1 + tmp2;  
                                         EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=-2*A_00*w30 - 2*A_11*w31 + A_22*w32 + tmp13 + tmp17 + tmp5;  
                                         EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=2*A_00*w30 + 32*A_11*w22 + 2*A_22*w32 + tmp11 + tmp15 + tmp7;  
                                         EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=-4*A_00*w30 + 2*A_11*w31 + 2*A_22*w32 + tmp16 + tmp8 + tmp9;  
                                         EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=4*A_00*w30 + 4*A_11*w31 + 4*A_22*w32 + tmp14 + tmp21 + tmp23;  
6076                                      }                                      }
6077                                  }                                  }
6078                              }                              }
# Line 6471  void Brick::assemblePDESystem(Paso_Syste Line 6443  void Brick::assemblePDESystem(Paso_Syste
6443                              } else { // constant data                              } else { // constant data
6444                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
6445                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
6446                                          const double B_0 = B_p[INDEX3(k,0,m,numEq,3)];                                          const double wB0 = B_p[INDEX3(k,0,m,numEq,3)]*w55;
6447                                          const double B_1 = B_p[INDEX3(k,1,m,numEq,3)];                                          const double wB1 = B_p[INDEX3(k,1,m,numEq,3)]*w56;
6448                                          const double B_2 = B_p[INDEX3(k,2,m,numEq,3)];                                          const double wB2 = B_p[INDEX3(k,2,m,numEq,3)]*w54;
6449                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=4*B_0*w55 + 4*B_1*w56 - B_2*(-4*w54);                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+= 4*wB0 + 4*wB1 + 4*wB2;
6450                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=4*B_0*w55 - B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+= 4*wB0 + 2*wB1 + 2*wB2;
6451                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=-B_0*(-2*w55) + 4*B_1*w56 + 2*B_2*w54;                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+= 2*wB0 + 4*wB1 + 2*wB2;
6452                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=-B_0*(-2*w55) - B_1*(-2*w56) + B_2*w54;                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+= 2*wB0 + 2*wB1 +   wB2;
6453                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=-B_0*(-2*w55) - B_1*(-2*w56) - B_2*(-4*w54);                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+= 2*wB0 + 2*wB1 + 4*wB2;
6454                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=-B_0*(-2*w55) + B_1*w56 + 2*B_2*w54;                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+= 2*wB0 +   wB1 + 2*wB2;
6455                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=B_0*w55 - B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=   wB0 + 2*wB1 + 2*wB2;
6456                                          EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=B_0*w55 + B_1*w56 + B_2*w54;                                          EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=   wB0 +   wB1 +   wB2;
6457                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=2*B_0*(-2*w55) - B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=-4*wB0 + 2*wB1 + 2*wB2;
6458                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=2*B_0*(-2*w55) + 4*B_1*w56 - B_2*(-4*w54);                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=-4*wB0 + 4*wB1 + 4*wB2;
6459                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=B_0*(-2*w55) - B_1*(-2*w56) + B_2*w54;                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=-2*wB0 + 2*wB1 +   wB2;
6460                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=B_0*(-2*w55) - 2*B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=-2*wB0 + 4*wB1 + 2*wB2;
6461                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=B_0*(-2*w55) + B_1*w56 + 2*B_2*w54;                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=-2*wB0 +   wB1 + 2*wB2;
6462                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=B_0*(-2*w55) - B_1*(-2*w56) - B_2*(-4*w54);                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=-2*wB0 + 2*wB1 + 4*wB2;
6463                                          EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=-B_0*w55 + B_1*w56 + B_2*w54;                                          EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=  -wB0 +   wB1 +   wB2;
6464                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=-B_0*w55 - B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=  -wB0 + 2*wB1 + 2*wB2;
6465                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=-B_0*(-2*w55) + 2*B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+= 2*wB0 - 4*wB1 + 2*wB2;
6466                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=-B_0*(-2*w55) + B_1*(-2*w56) + B_2*w54;                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+= 2*wB0 - 2*wB1 +   wB2;
6467                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=4*B_0*w55 + 2*B_1*(-2*w56) - B_2*(-4*w54);                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+= 4*wB0 - 4*wB1 + 4*wB2;
6468                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=4*B_0*w55 + B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+= 4*wB0 - 2*wB1 + 2*wB2;
6469                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=B_0*w55 + B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=   wB0 - 2*wB1 + 2*wB2;
6470                                          EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=B_0*w55 - B_1*w56 + B_2*w54;                                          EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=   wB0 -   wB1 +   wB2;
6471                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=-B_0*(-2*w55) + B_1*(-2*w56) - B_2*(-4*w54);                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+= 2*wB0 - 2*wB1 + 4*wB2;
6472                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=-B_0*(-2*w55) - B_1*w56 + 2*B_2*w54;                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+= 2*wB0 -   wB1 + 2*wB2;
6473                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=B_0*(-2*w55) + B_1*(-2*w56) + B_2*w54;                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=-2*wB0 - 2*wB1 +   wB2;
6474                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=B_0*(-2*w55) + 2*B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=-2*wB0 - 4*wB1 + 2*wB2;
6475                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=2*B_0*(-2*w55) + B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=-4*wB0 - 2*wB1 + 2*wB2;
6476                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=2*B_0*(-2*w55) + 2*B_1*(-2*w56) - B_2*(-4*w54);                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=-4*wB0 - 4*wB1 + 4*wB2;
6477                                          EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=-B_0*w55 - B_1*w56 + B_2*w54;                                          EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=  -wB0 -   wB1 +   wB2;
6478                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=-B_0*w55 + B_1*(-2*w56) + 2*B_2*w54;                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=  -wB0 - 2*wB1 + 2*wB2;
6479                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=B_0*(-2*w55) - B_1*w56 + 2*B_2*w54;                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=-2*wB0 -   wB1 + 2*wB2;
6480                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=B_0*(-2*w55) + B_1*(-2*w56) - B_2*(-4*w54);                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=-2*wB0 - 2*wB1 + 4*wB2;
6481                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=-B_0*(-2*w55) - B_1*(-2*w56) + B_2*(-4*w54);                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+= 2*wB0 + 2*wB1 - 4*wB2;
6482                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=-B_0*(-2*w55) + B_1*w56 - 2*B_2*w54;                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+= 2*wB0 +   wB1 - 2*wB2;
6483                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=B_0*w55 - B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=   wB0 + 2*wB1 - 2*wB2;
6484                                          EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=B_0*w55 + B_1*w56 + B_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=   wB0 +   wB1 -   wB2;
6485                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=4*B_0*w55 + 4*B_1*w56 + B_2*(-4*w54);                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+= 4*wB0 + 4*wB1 - 4*wB2;
6486                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=4*B_0*w55 - B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+= 4*wB0 + 2*wB1 - 2*wB2;
6487                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=-B_0*(-2*w55) + 4*B_1*w56 - 2*B_2*w54;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+= 2*wB0 + 4*wB1 - 2*wB2;
6488                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=-B_0*(-2*w55) - B_1*(-2*w56) + B_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+= 2*wB0 + 2*wB1 -   wB2;
6489                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=B_0*(-2*w55) + B_1*w56 - 2*B_2*w54;                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=-2*wB0 +   wB1 - 2*wB2;
6490                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=B_0*(-2*w55) - B_1*(-2*w56) + B_2*(-4*w54);                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=-2*wB0 + 2*wB1 - 4*wB2;
6491                                          EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=-B_0*w55 + B_1*w56 + B_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=  -wB0 +   wB1 -   wB2;
6492                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=-B_0*w55 - B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=  -wB0 + 2*wB1 - 2*wB2;
6493                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=2*B_0*(-2*w55) - B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=-4*wB0 + 2*wB1 - 2*wB2;
6494                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=2*B_0*(-2*w55) + 4*B_1*w56 + B_2*(-4*w54);                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=-4*wB0 + 4*wB1 - 4*wB2;
6495                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=B_0*(-2*w55) - B_1*(-2*w56) + B_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=-2*wB0 + 2*wB1 -   wB2;
6496                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=B_0*(-2*w55) + 4*B_1*w56 - 2*B_2*w54;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=-2*wB0 + 4*wB1 - 2*wB2;
6497                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=B_0*w55 + B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=   wB0 - 2*wB1 - 2*wB2;
6498                                          EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=B_0*w55 - B_1*w56 + B_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=   wB0 -   wB1 -   wB2;
6499                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=-B_0*(-2*w55) + B_1*(-2*w56) + B_2*(-4*w54);                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+= 2*wB0 - 2*wB1 - 4*wB2;
6500                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=-B_0*(-2*w55) - B_1*w56 - 2*B_2*w54;                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+= 2*wB0 -   wB1 - 2*wB2;
6501                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=-B_0*(-2*w55) + 2*B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+= 2*wB0 - 4*wB1 - 2*wB2;
6502                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=-B_0*(-2*w55) + B_1*(-2*w56) + B_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+= 2*wB0 - 2*wB1 -   wB2;
6503                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=4*B_0*w55 + 2*B_1*(-2*w56) + B_2*(-4*w54);                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+= 4*wB0 - 4*wB1 - 4*wB2;
6504                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=4*B_0*w55 + B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+= 4*wB0 - 2*wB1 - 2*wB2;
6505                                          EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=-B_0*w55 - B_1*w56 + B_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=  -wB0 -   wB1 -   wB2;
6506                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=-B_0*w55 + B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=  -wB0 - 2*wB1 - 2*wB2;
6507                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=B_0*(-2*w55) - B_1*w56 - 2*B_2*w54;                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=-2*wB0 -   wB1 - 2*wB2;
6508                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=B_0*(-2*w55) + B_1*(-2*w56) + B_2*(-4*w54);                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=-2*wB0 - 2*wB1 - 4*wB2;
6509                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=B_0*(-2*w55) + B_1*(-2*w56) + B_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=-2*wB0 - 2*wB1 -   wB2;
6510                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=B_0*(-2*w55) + 2*B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=-2*wB0 - 4*wB1 - 2*wB2;
6511                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=2*B_0*(-2*w55) + B_1*(-2*w56) - 2*B_2*w54;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=-4*wB0 - 2*wB1 - 2*wB2;
6512                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=2*B_0*(-2*w55) + 2*B_1*(-2*w56) + B_2*(-4*w54);                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=-4*wB0 - 4*wB1 - 4*wB2;
6513                                      }                                      }
6514                                  }                                  }
6515                              }                              }
# Line 6908  void Brick::assemblePDESystem(Paso_Syste Line 6880  void Brick::assemblePDESystem(Paso_Syste
6880                              } else { // constant data                              } else { // constant data
6881                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
6882                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
6883                                          const double C_0 = C_p[INDEX3(k,m,0,numEq,numComp)];                                          const double wC0 = C_p[INDEX3(k,m,0,numEq,numComp)]*w55;
6884                                          const double C_1 = C_p[INDEX3(k,m,1,numEq,numComp)];                                          const double wC1 = C_p[INDEX3(k,m,1,numEq,numComp)]*w56;
6885                                          const double C_2 = C_p[INDEX3(k,m,2,numEq,numComp)];                                          const double wC2 = C_p[INDEX3(k,m,2,numEq,numComp)]*w54;
6886                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=4*C_0*w55 + 4*C_1*w56 - C_2*(-4*w54);                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+= 4*wC0 + 4*wC1 + 4*wC2;
6887                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=2*C_0*(-2*w55) - C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=-4*wC0 + 2*wC1 + 2*wC2;
6888                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=-C_0*(-2*w55) + 2*C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+= 2*wC0 - 4*wC1 + 2*wC2;
6889                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=C_0*(-2*w55) + C_1*(-2*w56) + C_2*w54;                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=-2*wC0 - 2*wC1 +   wC2;
6890                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=-C_0*(-2*w55) - C_1*(-2*w56) + C_2*(-4*w54);                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+= 2*wC0 + 2*wC1 - 4*wC2;
6891                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=C_0*(-2*w55) + C_1*w56 - 2*C_2*w54;                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=-2*wC0 +   wC1 - 2*wC2;
6892                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=C_0*w55 + C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=   wC0 - 2*wC1 - 2*wC2;
6893                                          EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=-C_0*w55 - C_1*w56 + C_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=  -wC0 -   wC1 -   wC2;
6894                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=4*C_0*w55 - C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+= 4*wC0 + 2*wC1 + 2*wC2;
6895                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=2*C_0*(-2*w55) + 4*C_1*w56 - C_2*(-4*w54);                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=-4*wC0 + 4*wC1 + 4*wC2;
6896                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=-C_0*(-2*w55) + C_1*(-2*w56) + C_2*w54;                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+= 2*wC0 - 2*wC1 +   wC2;
6897                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=C_0*(-2*w55) + 2*C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=-2*wC0 - 4*wC1 + 2*wC2;
6898                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=-C_0*(-2*w55) + C_1*w56 - 2*C_2*w54;                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+= 2*wC0 +   wC1 - 2*wC2;
6899                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=C_0*(-2*w55) - C_1*(-2*w56) + C_2*(-4*w54);                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=-2*wC0 + 2*wC1 - 4*wC2;
6900                                          EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=C_0*w55 - C_1*w56 + C_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=   wC0 -   wC1 -   wC2;
6901                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=-C_0*w55 + C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=  -wC0 - 2*wC1 - 2*wC2;
6902                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=-C_0*(-2*w55) + 4*C_1*w56 + 2*C_2*w54;                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+= 2*wC0 + 4*wC1 + 2*wC2;
6903                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=C_0*(-2*w55) - C_1*(-2*w56) + C_2*w54;                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=-2*wC0 + 2*wC1 +   wC2;
6904                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=4*C_0*w55 + 2*C_1*(-2*w56) - C_2*(-4*w54);                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+= 4*wC0 - 4*wC1 + 4*wC2;
6905                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=2*C_0*(-2*w55) + C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=-4*wC0 - 2*wC1 + 2*wC2;
6906                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=C_0*w55 - C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=   wC0 + 2*wC1 - 2*wC2;
6907                                          EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=-C_0*w55 + C_1*w56 + C_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=  -wC0 +   wC1 -   wC2;
6908                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=-C_0*(-2*w55) + C_1*(-2*w56) + C_2*(-4*w54);                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+= 2*wC0 - 2*wC1 - 4*wC2;
6909                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=C_0*(-2*w55) - C_1*w56 - 2*C_2*w54;                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=-2*wC0 -   wC1 - 2*wC2;
6910                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=-C_0*(-2*w55) - C_1*(-2*w56) + C_2*w54;                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+= 2*wC0 + 2*wC1 +   wC2;
6911                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=C_0*(-2*w55) + 4*C_1*w56 + 2*C_2*w54;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=-2*wC0 + 4*wC1 + 2*wC2;
6912                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=4*C_0*w55 + C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+= 4*wC0 - 2*wC1 + 2*wC2;
6913                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=2*C_0*(-2*w55) + 2*C_1*(-2*w56) - C_2*(-4*w54);                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=-4*wC0 - 4*wC1 + 4*wC2;
6914                                          EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=C_0*w55 + C_1*w56 + C_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=   wC0 +   wC1 -   wC2;
6915                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=-C_0*w55 - C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=  -wC0 + 2*wC1 - 2*wC2;
6916                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=-C_0*(-2*w55) - C_1*w56 - 2*C_2*w54;                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+= 2*wC0 -   wC1 - 2*wC2;
6917                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=C_0*(-2*w55) + C_1*(-2*w56) + C_2*(-4*w54);                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=-2*wC0 - 2*wC1 - 4*wC2;
6918                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=-C_0*(-2*w55) - C_1*(-2*w56) - C_2*(-4*w54);                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+= 2*wC0 + 2*wC1 + 4*wC2;
6919                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=C_0*(-2*w55) + C_1*w56 + 2*C_2*w54;                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=-2*wC0 +   wC1 + 2*wC2;
6920                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=C_0*w55 + C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=   wC0 - 2*wC1 + 2*wC2;
6921                                          EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=-C_0*w55 - C_1*w56 + C_2*w54;                                          EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=  -wC0 -   wC1 +   wC2;
6922                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=4*C_0*w55 + 4*C_1*w56 + C_2*(-4*w54);                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+= 4*wC0 + 4*wC1 - 4*wC2;
6923                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=2*C_0*(-2*w55) - C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=-4*wC0 + 2*wC1 - 2*wC2;
6924                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=-C_0*(-2*w55) + 2*C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+= 2*wC0 - 4*wC1 - 2*wC2;
6925                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=C_0*(-2*w55) + C_1*(-2*w56) + C_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=-2*wC0 - 2*wC1 -   wC2;
6926                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=-C_0*(-2*w55) + C_1*w56 + 2*C_2*w54;                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+= 2*wC0 +   wC1 + 2*wC2;
6927                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=C_0*(-2*w55) - C_1*(-2*w56) - C_2*(-4*w54);                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=-2*wC0 + 2*wC1 + 4*wC2;
6928                                          EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=C_0*w55 - C_1*w56 + C_2*w54;                                          EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=   wC0 -   wC1 +   wC2;
6929                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=-C_0*w55 + C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=  -wC0 - 2*wC1 + 2*wC2;
6930                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=4*C_0*w55 - C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+= 4*wC0 + 2*wC1 - 2*wC2;
6931                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=2*C_0*(-2*w55) + 4*C_1*w56 + C_2*(-4*w54);                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=-4*wC0 + 4*wC1 - 4*wC2;
6932                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=-C_0*(-2*w55) + C_1*(-2*w56) + C_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+= 2*wC0 - 2*wC1 -   wC2;
6933                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=C_0*(-2*w55) + 2*C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=-2*wC0 - 4*wC1 - 2*wC2;
6934                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=C_0*w55 - C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=   wC0 + 2*wC1 + 2*wC2;
6935                                          EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=-C_0*w55 + C_1*w56 + C_2*w54;                                          EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=  -wC0 +   wC1 +   wC2;
6936                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=-C_0*(-2*w55) + C_1*(-2*w56) - C_2*(-4*w54);                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+= 2*wC0 - 2*wC1 + 4*wC2;
6937                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=C_0*(-2*w55) - C_1*w56 + 2*C_2*w54;                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=-2*wC0 -   wC1 + 2*wC2;
6938                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=-C_0*(-2*w55) + 4*C_1*w56 - 2*C_2*w54;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+= 2*wC0 + 4*wC1 - 2*wC2;
6939                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=C_0*(-2*w55) - C_1*(-2*w56) + C_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=-2*wC0 + 2*wC1 -   wC2;
6940                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=4*C_0*w55 + 2*C_1*(-2*w56) + C_2*(-4*w54);                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+= 4*wC0 - 4*wC1 - 4*wC2;
6941                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=2*C_0*(-2*w55) + C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=-4*wC0 - 2*wC1 - 2*wC2;
6942                                          EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=C_0*w55 + C_1*w56 + C_2*w54;                                          EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=   wC0 +   wC1 +   wC2;
6943                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=-C_0*w55 - C_1*(-2*w56) + 2*C_2*w54;                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=  -wC0 + 2*wC1 + 2*wC2;
6944                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=-C_0*(-2*w55) - C_1*w56 + 2*C_2*w54;                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+= 2*wC0 -   wC1 + 2*wC2;
6945                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=C_0*(-2*w55) + C_1*(-2*w56) - C_2*(-4*w54);                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=-2*wC0 - 2*wC1 + 4*wC2;
6946                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=-C_0*(-2*w55) - C_1*(-2*w56) + C_2*(-4*w54)/4;                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+= 2*wC0 + 2*wC1 -   wC2;
6947                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=C_0*(-2*w55) + 4*C_1*w56 - 2*C_2*w54;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=-2*wC0 + 4*wC1 - 2*wC2;
6948                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=4*C_0*w55 + C_1*(-2*w56) - 2*C_2*w54;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+= 4*wC0 - 2*wC1 - 2*wC2;
6949                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=2*C_0*(-2*w55) + 2*C_1*(-2*w56) + C_2*(-4*w54);                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=-4*wC0 - 4*wC1 - 4*wC2;
6950                                      }                                      }
6951                                  }                                  }
6952                              }                              }
# Line 7124  void Brick::assemblePDESystem(Paso_Syste Line 7096  void Brick::assemblePDESystem(Paso_Syste
7096                              } else { // constant data                              } else { // constant data
7097                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
7098                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
7099                                          const double D_0 = D_p[INDEX2(k, m, numEq)];                                          const double wD0 = 8*D_p[INDEX2(k, m, numEq)]*w58;
7100                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=64*D_0*w58;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=8*wD0;
7101                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=4*wD0;
7102                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=4*wD0;
7103                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=2*wD0;
7104                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=4*wD0;
7105                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=2*wD0;
7106                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=2*wD0;
7107                                          EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=8*D_0*w58;                                          EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=  wD0;
7108                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=4*wD0;
7109                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=64*D_0*w58;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=8*wD0;
7110                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=2*wD0;
7111                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=4*wD0;
7112                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=2*wD0;
7113                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=4*wD0;
7114                                          EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=8*D_0*w58;                                          EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=  wD0;
7115                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=2*wD0;
7116                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=4*wD0;
7117                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=2*wD0;
7118                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=64*D_0*w58;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=8*wD0;
7119                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=4*wD0;
7120                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=2*wD0;
7121                                          EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=8*D_0*w58;                                          EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=  wD0;
7122                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=4*wD0;
7123                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=2*wD0;
7124                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=2*wD0;
7125                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=4*wD0;
7126                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=4*wD0;
7127                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=64*D_0*w58;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=8*wD0;
7128                                          EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=8*D_0*w58;                                          EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=  wD0;
7129                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=2*wD0;
7130                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=2*wD0;
7131                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=4*wD0;
7132                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=4*wD0;
7133                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=2*wD0;
7134                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=2*wD0;
7135                                          EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=8*D_0*w58;                                          EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=  wD0;
7136                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=64*D_0*w58;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=8*wD0;
7137                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=4*wD0;
7138                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=4*wD0;
7139                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=2*wD0;
7140                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=2*wD0;
7141                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=4*wD0;
7142                                          EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=8*D_0*w58;                                          EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=  wD0;
7143                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=2*wD0;
7144                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=4*wD0;
7145                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=64*D_0*w58;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=8*wD0;
7146                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=2*wD0;
7147                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=4*wD0;
7148                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=2*wD0;
7149                                          EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=8*D_0*w58;                                          EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=  wD0;
7150                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=4*wD0;
7151                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=2*wD0;
7152                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=4*wD0;
7153                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=2*wD0;
7154                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=64*D_0*w58;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=8*wD0;
7155                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=4*wD0;
7156                                          EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=8*D_0*w58;                                          EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=  wD0;
7157                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=2*wD0;
7158                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=2*wD0;
7159                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=4*wD0;
7160                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=16*D_0*w58;                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=2*wD0;
7161                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=4*wD0;
7162                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=32*D_0*w58;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=4*wD0;
7163                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=64*D_0*w58;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=8*wD0;
7164                                      }                                      }
7165                                  }                                  }
7166                              }                              }
# Line 7296  void Brick::assemblePDESystem(Paso_Syste Line 7268  void Brick::assemblePDESystem(Paso_Syste
7268                                  }                                  }
7269                              } else { // constant data                              } else { // constant data
7270                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
7271                                      const double X_0 = X_p[INDEX2(k, 0, numEq)];                                      const double wX0 = 18*X_p[INDEX2(k, 0, numEq)]*w55;
7272                                      const double X_1 = X_p[INDEX2(k, 1, numEq)];                                      const double wX1 = 18*X_p[INDEX2(k, 1, numEq)]*w56;
7273                                      const double X_2 = X_p[INDEX2(k, 2, numEq)];                                      const double wX2 = 18*X_p[INDEX2(k, 2, numEq)]*w54;
7274                                      EM_F[INDEX2(k,0,numEq)]+=12*X_0*w69 + 12*X_1*w67 + 18*X_2*w54;                                      EM_F[INDEX2(k,0,numEq)]+= wX0 + wX1 + wX2;
7275                                      EM_F[INDEX2(k,1,numEq)]+=-18*X_0*w55 + 12*X_1*w67 + 18*X_2*w54;                                      EM_F[INDEX2(k,1,numEq)]+=-wX0 + wX1 + wX2;
7276                                      EM_F[INDEX2(k,2,numEq)]+=12*X_0*w69 - 18*X_1*w56 + 18*X_2*w54;                                      EM_F[INDEX2(k,2,numEq)]+= wX0 - wX1 + wX2;
7277                                      EM_F[INDEX2(k,3,numEq)]+=-18*X_0*w55 - 18*X_1*w56 + 18*X_2*w54;                                      EM_F[INDEX2(k,3,numEq)]+=-wX0 - wX1 + wX2;
7278                                      EM_F[INDEX2(k,4,numEq)]+=12*X_0*w69 + 12*X_1*w67 - 18*X_2*w54;                                      EM_F[INDEX2(k,4,numEq)]+= wX0 + wX1 - wX2;
7279                                      EM_F[INDEX2(k,5,numEq)]+=-18*X_0*w55 + 12*X_1*w67 - 18*X_2*w54;                                      EM_F[INDEX2(k,5,numEq)]+=-wX0 + wX1 - wX2;
7280                                      EM_F[INDEX2(k,6,numEq)]+=12*X_0*w69 - 18*X_1*w56 - 18*X_2*w54;                                      EM_F[INDEX2(k,6,numEq)]+= wX0 - wX1 - wX2;
7281                                      EM_F[INDEX2(k,7,numEq)]+=-18*X_0*w55 - 18*X_1*w56 - 18*X_2*w54;                                      EM_F[INDEX2(k,7,numEq)]+=-wX0 - wX1 - wX2;
7282                                  }                                  }
7283                              }                              }
7284                          }                          }
# Line 7785  void Brick::assemblePDEBoundarySingle(Pa Line 7757  void Brick::assemblePDEBoundarySingle(Pa
7757  {  {
7758      const double SQRT3 = 1.73205080756887719318;      const double SQRT3 = 1.73205080756887719318;
7759      const double w12 = m_dx[0]*m_dx[1]/144;      const double w12 = m_dx[0]*m_dx[1]/144;
7760        const double w10 = w12*(-SQRT3 + 2);
7761      const double w11 = w12*(SQRT3 + 2);      const double w11 = w12*(SQRT3 + 2);
     const double w14 = w12*(4*SQRT3 + 7);  
7762      const double w13 = w12*(-4*SQRT3 + 7);      const double w13 = w12*(-4*SQRT3 + 7);
7763      const double w10 = w12*(-SQRT3 + 2);      const double w14 = w12*(4*SQRT3 + 7);
7764      const double w7 = m_dx[0]*m_dx[2]/144;      const double w7 = m_dx[0]*m_dx[2]/144;
7765        const double w5 = w7*(-SQRT3 + 2);
7766      const double w6 = w7*(SQRT3 + 2);      const double w6 = w7*(SQRT3 + 2);
7767      const double w8 = w7*(-4*SQRT3 + 7);      const double w8 = w7*(-4*SQRT3 + 7);
7768      const double w9 = w7*(4*SQRT3 + 7);      const double w9 = w7*(4*SQRT3 + 7);
     const double w5 = w7*(-SQRT3 + 2);  
7769      const double w2 = m_dx[1]*m_dx[2]/144;      const double w2 = m_dx[1]*m_dx[2]/144;
7770        const double w0 = w2*(-SQRT3 + 2);
7771      const double w1 = w2*(SQRT3 + 2);      const double w1 = w2*(SQRT3 + 2);
     const double w4 = w2*(4*SQRT3 + 7);  
7772      const double w3 = w2*(-4*SQRT3 + 7);      const double w3 = w2*(-4*SQRT3 + 7);
7773      const double w0 = w2*(-SQRT3 + 2);      const double w4 = w2*(4*SQRT3 + 7);
7774      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
7775      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
7776      rhs.requireWrite();      rhs.requireWrite();
# Line 7850  void Brick::assemblePDEBoundarySingle(Pa Line 7822  void Brick::assemblePDEBoundarySingle(Pa
7822                                  EM_S[INDEX2(6,4,8)]+=tmp0 + tmp1;                                  EM_S[INDEX2(6,4,8)]+=tmp0 + tmp1;
7823                                  EM_S[INDEX2(6,6,8)]+=d_0*w3 + d_3*w4 + tmp9;                                  EM_S[INDEX2(6,6,8)]+=d_0*w3 + d_3*w4 + tmp9;
7824                              } else { // constant data                              } else { // constant data
7825                                  const double d_0 = d_p[0];                                  const double wd0 = 4*d_p[0]*w2;
7826                                  EM_S[INDEX2(0,0,8)]+=16*d_0*w2;                                  EM_S[INDEX2(0,0,8)]+=4*wd0;
7827                                  EM_S[INDEX2(0,2,8)]+=8*d_0*w2;                                  EM_S[INDEX2(0,2,8)]+=2*wd0;
7828                                  EM_S[INDEX2(0,4,8)]+=8*d_0*w2;                                  EM_S[INDEX2(0,4,8)]+=2*wd0;
7829                                  EM_S[INDEX2(0,6,8)]+=4*d_0*w2;                                  EM_S[INDEX2(0,6,8)]+=  wd0;
7830                                  EM_S[INDEX2(2,0,8)]+=8*d_0*w2;                                  EM_S[INDEX2(2,0,8)]+=2*wd0;
7831                                  EM_S[INDEX2(2,2,8)]+=16*d_0*w2;                                  EM_S[INDEX2(2,2,8)]+=4*wd0;
7832                                  EM_S[INDEX2(2,4,8)]+=4*d_0*w2;                                  EM_S[INDEX2(2,4,8)]+=  wd0;
7833                                  EM_S[INDEX2(2,6,8)]+=8*d_0*w2;                                  EM_S[INDEX2(2,6,8)]+=2*wd0;
7834                                  EM_S[INDEX2(4,0,8)]+=8*d_0*w2;                                  EM_S[INDEX2(4,0,8)]+=2*wd0;
7835                                  EM_S[INDEX2(4,2,8)]+=4*d_0*w2;                                  EM_S[INDEX2(4,2,8)]+=  wd0;
7836                                  EM_S[INDEX2(4,4,8)]+=16*d_0*w2;                                  EM_S[INDEX2(4,4,8)]+=4*wd0;
7837                                  EM_S[INDEX2(4,6,8)]+=8*d_0*w2;                                  EM_S[INDEX2(4,6,8)]+=2*wd0;
7838                                  EM_S[INDEX2(6,0,8)]+=4*d_0*w2;                                  EM_S[INDEX2(6,0,8)]+=  wd0;
7839                                  EM_S[INDEX2(6,2,8)]+=8*d_0*w2;                                  EM_S[INDEX2(6,2,8)]+=2*wd0;
7840                                  EM_S[INDEX2(6,4,8)]+=8*d_0*w2;                                  EM_S[INDEX2(6,4,8)]+=2*wd0;
7841                                  EM_S[INDEX2(6,6,8)]+=16*d_0*w2;                                  EM_S[INDEX2(6,6,8)]+=4*wd0;
7842                              }                              }
7843                          }                          }
7844                          ///////////////                          ///////////////
# Line 7946  void Brick::assemblePDEBoundarySingle(Pa Line 7918  void Brick::assemblePDEBoundarySingle(Pa
7918                                  EM_S[INDEX2(7,5,8)]+=tmp8 + tmp9;                                  EM_S[INDEX2(7,5,8)]+=tmp8 + tmp9;
7919                                  EM_S[INDEX2(7,7,8)]+=d_0*w3 + d_3*w4 + tmp7;                                  EM_S[INDEX2(7,7,8)]+=d_0*w3 + d_3*w4 + tmp7;
7920                              } else { // constant data                              } else { // constant data
7921                                  const double d_0 = d_p[0];                                  const double wd0 = 4*d_p[0]*w2;
7922                                  EM_S[INDEX2(1,1,8)]+=16*d_0*w2;                                  EM_S[INDEX2(1,1,8)]+=4*wd0;
7923                                  EM_S[INDEX2(1,3,8)]+=8*d_0*w2;                                  EM_S[INDEX2(1,3,8)]+=2*wd0;
7924                                  EM_S[INDEX2(1,5,8)]+=8*d_0*w2;                                  EM_S[INDEX2(1,5,8)]+=2*wd0;
7925                                  EM_S[INDEX2(1,7,8)]+=4*d_0*w2;                                  EM_S[INDEX2(1,7,8)]+=  wd0;
7926                                  EM_S[INDEX2(3,1,8)]+=8*d_0*w2;                                  EM_S[INDEX2(3,1,8)]+=2*wd0;
7927                                  EM_S[INDEX2(3,3,8)]+=16*d_0*w2;                                  EM_S[INDEX2(3,3,8)]+=4*wd0;
7928                                  EM_S[INDEX2(3,5,8)]+=4*d_0*w2;                                  EM_S[INDEX2(3,5,8)]+=  wd0;
7929                                  EM_S[INDEX2(3,7,8)]+=8*d_0*w2;                                  EM_S[INDEX2(3,7,8)]+=2*wd0;
7930                                  EM_S[INDEX2(5,1,8)]+=8*d_0*w2;                                  EM_S[INDEX2(5,1,8)]+=2*wd0;
7931                                  EM_S[INDEX2(5,3,8)]+=4*d_0*w2;                                  EM_S[INDEX2(5,3,8)]+=  wd0;
7932                                  EM_S[INDEX2(5,5,8)]+=16*d_0*w2;                                  EM_S[INDEX2(5,5,8)]+=4*wd0;
7933                                  EM_S[INDEX2(5,7,8)]+=8*d_0*w2;                                  EM_S[INDEX2(5,7,8)]+=2*wd0;
7934                                  EM_S[INDEX2(7,1,8)]+=4*d_0*w2;                                  EM_S[INDEX2(7,1,8)]+=  wd0;
7935                                  EM_S[INDEX2(7,3,8)]+=8*d_0*w2;                                  EM_S[INDEX2(7,3,8)]+=2*wd0;
7936                                  EM_S[INDEX2(7,5,8)]+=8*d_0*w2;                                  EM_S[INDEX2(7,5,8)]+=2*wd0;
7937                                  EM_S[INDEX2(7,7,8)]+=16*d_0*w2;                                  EM_S[INDEX2(7,7,8)]+=4*wd0;
7938                              }                              }
7939                          }                          }
7940                          ///////////////                          ///////////////
# Line 8042  void Brick::assemblePDEBoundarySingle(Pa Line 8014  void Brick::assemblePDEBoundarySingle(Pa
8014                                  EM_S[INDEX2(5,4,8)]+=tmp0 + tmp1;                                  EM_S[INDEX2(5,4,8)]+=tmp0 + tmp1;
8015                                  EM_S[INDEX2(5,5,8)]+=d_0*w8 + d_3*w9 + tmp8;                                  EM_S[INDEX2(5,5,8)]+=d_0*w8 + d_3*w9 + tmp8;
8016                              } else { // constant data                              } else { // constant data
8017                                  const double d_0 = d_p[0];                                  const double wd0 = 4*d_p[0]*w7;
8018                                  EM_S[INDEX2(0,0,8)]+=16*d_0*w7;                                  EM_S[INDEX2(0,0,8)]+=4*wd0;
8019                                  EM_S[INDEX2(0,1,8)]+=8*d_0*w7;                                  EM_S[INDEX2(0,1,8)]+=2*wd0;
8020                                  EM_S[INDEX2(0,4,8)]+=8*d_0*w7;                                  EM_S[INDEX2(0,4,8)]+=2*wd0;
8021                                  EM_S[INDEX2(0,5,8)]+=4*d_0*w7;                                  EM_S[INDEX2(0,5,8)]+=  wd0;
8022                                  EM_S[INDEX2(1,0,8)]+=8*d_0*w7;                                  EM_S[INDEX2(1,0,8)]+=2*wd0;
8023                                  EM_S[INDEX2(1,1,8)]+=16*d_0*w7;                                  EM_S[INDEX2(1,1,8)]+=4*wd0;
8024                                  EM_S[INDEX2(1,4,8)]+=4*d_0*w7;                                  EM_S[INDEX2(1,4,8)]+=  wd0;
8025                                  EM_S[INDEX2(1,5,8)]+=8*d_0*w7;                                  EM_S[INDEX2(1,5,8)]+=2*wd0;
8026                                  EM_S[INDEX2(4,0,8)]+=8*d_0*w7;                                  EM_S[INDEX2(4,0,8)]+=2*wd0;
8027                                  EM_S[INDEX2(4,1,8)]+=4*d_0*w7;                                  EM_S[INDEX2(4,1,8)]+=  wd0;
8028                                  EM_S[INDEX2(4,4,8)]+=16*d_0*w7;                                  EM_S[INDEX2(4,4,8)]+=4*wd0;
8029                                  EM_S[INDEX2(4,5,8)]+=8*d_0*w7;                                  EM_S[INDEX2(4,5,8)]+=2*wd0;
8030                                  EM_S[INDEX2(5,0,8)]+=4*d_0*w7;                                  EM_S[INDEX2(5,0,8)]+=  wd0;
8031                                  EM_S[INDEX2(5,1,8)]+=8*d_0*w7;                                  EM_S[INDEX2(5,1,8)]+=2*wd0;
8032                                  EM_S[INDEX2(5,4,8)]+=8*d_0*w7;                                  EM_S[INDEX2(5,4,8)]+=2*wd0;
8033                                  EM_S[INDEX2(5,5,8)]+=16*d_0*w7;                                  EM_S[INDEX2(5,5,8)]+=4*wd0;
8034                              }                              }
8035                          }                          }
8036                          ///////////////                          ///////////////
# Line 8138  void Brick::assemblePDEBoundarySingle(Pa Line 8110  void Brick::assemblePDEBoundarySingle(Pa
8110                                  EM_S[INDEX2(7,6,8)]+=tmp5 + tmp6;                                  EM_S[INDEX2(7,6,8)]+=tmp5 + tmp6;
8111                                  EM_S[INDEX2(7,7,8)]+=d_0*w8 + d_3*w9 + tmp8;                                  EM_S[INDEX2(7,7,8)]+=d_0*w8 + d_3*w9 + tmp8;
8112                              } else { // constant data                              } else { // constant data
8113                                  const double d_0 = d_p[0];                                  const double wd0 = 4*d_p[0]*w7;
8114                                  EM_S[INDEX2(2,2,8)]+=16*d_0*w7;                                  EM_S[INDEX2(2,2,8)]+=4*wd0;
8115                                  EM_S[INDEX2(2,3,8)]+=8*d_0*w7;                                  EM_S[INDEX2(2,3,8)]+=2*wd0;
8116                                  EM_S[INDEX2(2,6,8)]+=8*d_0*w7;                                  EM_S[INDEX2(2,6,8)]+=2*wd0;
8117                                  EM_S[INDEX2(2,7,8)]+=4*d_0*w7;                                  EM_S[INDEX2(2,7,8)]+=  wd0;
8118                                  EM_S[INDEX2(3,2,8)]+=8*d_0*w7;                                  EM_S[INDEX2(3,2,8)]+=2*wd0;
8119                                  EM_S[INDEX2(3,3,8)]+=16*d_0*w7;                                  EM_S[INDEX2(3,3,8)]+=4*wd0;
8120                                  EM_S[INDEX2(3,6,8)]+=4*d_0*w7;                                  EM_S[INDEX2(3,6,8)]+=  wd0;
8121                                  EM_S[INDEX2(3,7,8)]+=8*d_0*w7;                                  EM_S[INDEX2(3,7,8)]+=2*wd0;
8122                                  EM_S[INDEX2(6,2,8)]+=8*d_0*w7;                                  EM_S[INDEX2(6,2,8)]+=2*wd0;
8123                                  EM_S[INDEX2(6,3,8)]+=4*d_0*w7;                                  EM_S[INDEX2(6,3,8)]+=  wd0;
8124                                  EM_S[INDEX2(6,6,8)]+=16*d_0*w7;                                  EM_S[INDEX2(6,6,8)]+=4*wd0;
8125                                  EM_S[INDEX2(6,7,8)]+=8*d_0*w7;                                  EM_S[INDEX2(6,7,8)]+=2*wd0;
8126                                  EM_S[INDEX2(7,2,8)]+=4*d_0*w7;                                  EM_S[INDEX2(7,2,8)]+=  wd0;
8127                                  EM_S[INDEX2(7,3,8)]+=8*d_0*w7;                                  EM_S[INDEX2(7,3,8)]+=2*wd0;
8128                                  EM_S[INDEX2(7,6,8)]+=8*d_0*w7;                                  EM_S[INDEX2(7,6,8)]+=2*wd0;
8129                                  EM_S[INDEX2(7,7,8)]+=16*d_0*w7;                                  EM_S[INDEX2(7,7,8)]+=4*wd0;
8130                              }                              }
8131                          }                          }
8132                          ///////////////                          ///////////////
# Line 8234  void Brick::assemblePDEBoundarySingle(Pa Line 8206  void Brick::assemblePDEBoundarySingle(Pa
8206                                  EM_S[INDEX2(3,2,8)]+=tmp7 + tmp8;                                  EM_S[INDEX2(3,2,8)]+=tmp7 + tmp8;
8207                                  EM_S[INDEX2(3,3,8)]+=d_0*w13 + d_3*w14 + tmp3;                                  EM_S[INDEX2(3,3,8)]+=d_0*w13 + d_3*w14 + tmp3;
8208                              } else { // constant data                              } else { // constant data
8209                                  const double d_0 = d_p[0];                                  const double wd0 = 4*d_p[0]*w12;
8210                                  EM_S[INDEX2(0,0,8)]+=16*d_0*w12;                                  EM_S[INDEX2(0,0,8)]+=4*wd0;
8211                                  EM_S[INDEX2(0,1,8)]+=8*d_0*w12;                                  EM_S[INDEX2(0,1,8)]+=2*wd0;
8212                                  EM_S[INDEX2(0,2,8)]+=8*d_0*w12;                                  EM_S[INDEX2(0,2,8)]+=2*wd0;
8213                                  EM_S[INDEX2(0,3,8)]+=4*d_0*w12;                                  EM_S[INDEX2(0,3,8)]+=  wd0;
8214                                  EM_S[INDEX2(1,0,8)]+=8*d_0*w12;                                  EM_S[INDEX2(1,0,8)]+=2*wd0;
8215                                  EM_S[INDEX2(1,1,8)]+=16*d_0*w12;                                  EM_S[INDEX2(1,1,8)]+=4*wd0;
8216                                  EM_S[INDEX2(1,2,8)]+=4*d_0*w12;                                  EM_S[INDEX2(1,2,8)]+=  wd0;
8217                                  EM_S[INDEX2(1,3,8)]+=8*d_0*w12;                                  EM_S[INDEX2(1,3,8)]+=2*wd0;
8218                                  EM_S[INDEX2(2,0,8)]+=8*d_0*w12;                                  EM_S[INDEX2(2,0,8)]+=2*wd0;
8219                                  EM_S[INDEX2(2,1,8)]+=4*d_0*w12;                                  EM_S[INDEX2(2,1,8)]+=  wd0;
8220                                  EM_S[INDEX2(2,2,8)]+=16*d_0*w12;                                  EM_S[INDEX2(2,2,8)]+=4*wd0;
8221                                  EM_S[INDEX2(2,3,8)]+=8*d_0*w12;                                  EM_S[INDEX2(2,3,8)]+=2*wd0;
8222                                  EM_S[INDEX2(3,0,8)]+=4*d_0*w12;                                  EM_S[INDEX2(3,0,8)]+=  wd0;
8223                                  EM_S[INDEX2(3,1,8)]+=8*d_0*w12;                                  EM_S[INDEX2(3,1,8)]+=2*wd0;
8224                                  EM_S[INDEX2(3,2,8)]+=8*d_0*w12;                                  EM_S[INDEX2(3,2,8)]+=2*wd0;
8225                                  EM_S[INDEX2(3,3,8)]+=16*d_0*w12;                                  EM_S[INDEX2(3,3,8)]+=4*wd0;
8226                              }                              }
8227                          }                          }
8228                          ///////////////                          ///////////////
# Line 8330  void Brick::assemblePDEBoundarySingle(Pa Line 8302  void Brick::assemblePDEBoundarySingle(Pa
8302                                  EM_S[INDEX2(7,6,8)]+=tmp5 + tmp6;                                  EM_S[INDEX2(7,6,8)]+=tmp5 + tmp6;
8303                                  EM_S[INDEX2(7,7,8)]+=d_0*w13 + d_3*w14 + tmp7;                                  EM_S[INDEX2(7,7,8)]+=d_0*w13 + d_3*w14 + tmp7;
8304                              } else { // constant data                              } else { // constant data
8305                                  const double d_0 = d_p[0];                                  const double wd0 = 4*d_p[0]*w12;
8306                                  EM_S[INDEX2(4,4,8)]+=16*d_0*w12;                                  EM_S[INDEX2(4,4,8)]+=4*wd0;
8307                                  EM_S[INDEX2(4,5,8)]+=8*d_0*w12;                                  EM_S[INDEX2(4,5,8)]+=2*wd0;
8308                                  EM_S[INDEX2(4,6,8)]+=8*d_0*w12;                                  EM_S[INDEX2(4,6,8)]+=2*wd0;
8309                                  EM_S[INDEX2(4,7,8)]+=4*d_0*w12;                                  EM_S[INDEX2(4,7,8)]+=  wd0;
8310                                  EM_S[INDEX2(5,4,8)]+=8*d_0*w12;                                  EM_S[INDEX2(5,4,8)]+=2*wd0;
8311                                  EM_S[INDEX2(5,5,8)]+=16*d_0*w12;                                  EM_S[INDEX2(5,5,8)]+=4*wd0;
8312                                  EM_S[INDEX2(5,6,8)]+=4*d_0*w12;                                  EM_S[INDEX2(5,6,8)]+=  wd0;
8313                                  EM_S[INDEX2(5,7,8)]+=8*d_0*w12;                                  EM_S[INDEX2(5,7,8)]+=2*wd0;
8314                                  EM_S[INDEX2(6,4,8)]+=8*d_0*w12;                                  EM_S[INDEX2(6,4,8)]+=2*wd0;
8315                                  EM_S[INDEX2(6,5,8)]+=4*d_0*w12;                                  EM_S[INDEX2(6,5,8)]+=  wd0;
8316                                  EM_S[INDEX2(6,6,8)]+=16*d_0*w12;                                  EM_S[INDEX2(6,6,8)]+=4*wd0;
8317                                  EM_S[INDEX2(6,7,8)]+=8*d_0*w12;                                  EM_S[INDEX2(6,7,8)]+=2*wd0;
8318                                  EM_S[INDEX2(7,4,8)]+=4*d_0*w12;                                  EM_S[INDEX2(7,4,8)]+=  wd0;
8319                                  EM_S[INDEX2(7,5,8)]+=8*d_0*w12;                                  EM_S[INDEX2(7,5,8)]+=2*wd0;
8320                                  EM_S[INDEX2(7,6,8)]+=8*d_0*w12;                                  EM_S[INDEX2(7,6,8)]+=2*wd0;
8321                                  EM_S[INDEX2(7,7,8)]+=16*d_0*w12;                                  EM_S[INDEX2(7,7,8)]+=4*wd0;
8322                              }                              }
8323                          }                          }
8324                          ///////////////                          ///////////////
# Line 8386  void Brick::assemblePDEBoundarySingle(Pa Line 8358  void Brick::assemblePDEBoundarySingle(Pa
8358  void Brick::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,  void Brick::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
8359        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
8360  {  {
8361      const double w0 = m_dx[0]*m_dx[1]/16.;      const double w0 = m_dx[0]*m_dx[1]/16;
8362      const double w1 = m_dx[0]*m_dx[2]/16.;      const double w1 = m_dx[0]*m_dx[2]/16;
8363      const double w2 = m_dx[1]*m_dx[2]/16.;      const double w2 = m_dx[1]*m_dx[2]/16;
8364      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
8365      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
8366      rhs.requireWrite();      rhs.requireWrite();
# Line 8695  void Brick::assemblePDEBoundarySystem(Pa Line 8667  void Brick::assemblePDEBoundarySystem(Pa
8667          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
8668          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
8669      }      }
     /* GENERATOR SNIP_PDEBC_SYSTEM_PRE TOP */  
8670      const double SQRT3 = 1.73205080756887719318;      const double SQRT3 = 1.73205080756887719318;
8671      const double w12 = m_dx[0]*m_dx[1]/144;      const double w12 = m_dx[0]*m_dx[1]/144;
8672        const double w10 = w12*(-SQRT3 + 2);
8673      const double w11 = w12*(SQRT3 + 2);      const double w11 = w12*(SQRT3 + 2);
     const double w14 = w12*(4*SQRT3 + 7);  
8674      const double w13 = w12*(-4*SQRT3 + 7);      const double w13 = w12*(-4*SQRT3 + 7);
8675      const double w10 = w12*(-SQRT3 + 2);      const double w14 = w12*(4*SQRT3 + 7);
8676      const double w7 = m_dx[0]*m_dx[2]/144;      const double w7 = m_dx[0]*m_dx[2]/144;
8677        const double w5 = w7*(-SQRT3 + 2);
8678      const double w6 = w7*(SQRT3 + 2);      const double w6 = w7*(SQRT3 + 2);
8679      const double w8 = w7*(-4*SQRT3 + 7);      const double w8 = w7*(-4*SQRT3 + 7);
8680      const double w9 = w7*(4*SQRT3 + 7);      const double w9 = w7*(4*SQRT3 + 7);
     const double w5 = w7*(-SQRT3 + 2);  
8681      const double w2 = m_dx[1]*m_dx[2]/144;      const double w2 = m_dx[1]*m_dx[2]/144;
8682        const double w0 = w2*(-SQRT3 + 2);
8683      const double w1 = w2*(SQRT3 + 2);      const double w1 = w2*(SQRT3 + 2);
     const double w4 = w2*(4*SQRT3 + 7);  
8684      const double w3 = w2*(-4*SQRT3 + 7);      const double w3 = w2*(-4*SQRT3 + 7);
8685      const double w0 = w2*(-SQRT3 + 2);      const double w4 = w2*(4*SQRT3 + 7);
     /* GENERATOR SNIP_PDEBC_SYSTEM_PRE BOTTOM */  
8686      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
8687      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
8688      rhs.requireWrite();      rhs.requireWrite();
# Line 8726  void Brick::assemblePDEBoundarySystem(Pa Line 8696  void Brick::assemblePDEBoundarySystem(Pa
8696                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
8697                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
8698                          const index_t e = INDEX2(k1,k2,m_NE[1]);                          const index_t e = INDEX2(k1,k2,m_NE[1]);
                         /* GENERATOR SNIP_PDEBC_SYSTEM_0 TOP */  
8699                          ///////////////                          ///////////////
8700                          // process d //                          // process d //
8701                          ///////////////                          ///////////////
# Line 8771  void Brick::assemblePDEBoundarySystem(Pa Line 8740  void Brick::assemblePDEBoundarySystem(Pa
8740                              } else { // constant data                              } else { // constant data
8741                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8742                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8743                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double wd0 = 4*d_p[INDEX2(k, m, numEq)]*w2;
8744                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=16*d_0*w2;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=4*wd0;
8745                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=2*wd0;
8746                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=2*wd0;
8747                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=4*d_0*w2;                                          EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=  wd0;
8748                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=2*wd0;
8749                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=16*d_0*w2;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=4*wd0;
8750                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=4*d_0*w2;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=  wd0;
8751                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=2*wd0;
8752                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=2*wd0;
8753                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=4*d_0*w2;                                          EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=  wd0;
8754                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=16*d_0*w2;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=4*wd0;
8755                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=2*wd0;
8756                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=4*d_0*w2;                                          EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=  wd0;
8757                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=2*wd0;
8758                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=2*wd0;
8759                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=16*d_0*w2;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=4*wd0;
8760                                      }                                      }
8761                                  }                                  }
8762                              }                              }
# Line 8819  void Brick::assemblePDEBoundarySystem(Pa Line 8788  void Brick::assemblePDEBoundarySystem(Pa
8788                                  }                                  }
8789                              }                              }
8790                          }                          }
                         /* GENERATOR SNIP_PDEBC_SYSTEM_0 BOTTOM */  
8791                          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;
8792                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
8793                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 8836  void Brick::assemblePDEBoundarySystem(Pa Line 8804  void Brick::assemblePDEBoundarySystem(Pa
8804                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
8805                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
8806                          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]);
                         /* GENERATOR SNIP_PDEBC_SYSTEM_1 TOP */  
8807                          ///////////////                          ///////////////
8808                          // process d //                          // process d //
8809                          ///////////////                          ///////////////
# Line 8881  void Brick::assemblePDEBoundarySystem(Pa Line 8848  void Brick::assemblePDEBoundarySystem(Pa
8848                              } else { // constant data                              } else { // constant data
8849                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8850                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8851                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double wd0 = 4*d_p[INDEX2(k, m, numEq)]*w2;
8852                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=16*d_0*w2;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=4*wd0;
8853                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=2*wd0;
8854                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=2*wd0;
8855                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=4*d_0*w2;                                          EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=  wd0;
8856                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=2*wd0;
8857                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=16*d_0*w2;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=4*wd0;
8858                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=4*d_0*w2;                                          EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=  wd0;
8859                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=2*wd0;
8860                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=2*wd0;
8861                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=4*d_0*w2;                                          EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=  wd0;
8862                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=16*d_0*w2;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=4*wd0;
8863                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=2*wd0;
8864                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=4*d_0*w2;                                          EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=  wd0;
8865                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=2*wd0;
8866                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=8*d_0*w2;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=2*wd0;
8867                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=16*d_0*w2;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=4*wd0;
8868                                      }                                      }
8869                                  }                                  }
8870                              }                              }
# Line 8929  void Brick::assemblePDEBoundarySystem(Pa Line 8896  void Brick::assemblePDEBoundarySystem(Pa
8896                                  }                                  }
8897                              }                              }
8898                          }                          }
                         /* GENERATOR SNIP_PDEBC_SYSTEM_1 BOTTOM */  
8899                          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;
8900                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
8901                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 8946  void Brick::assemblePDEBoundarySystem(Pa Line 8912  void Brick::assemblePDEBoundarySystem(Pa
8912                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
8913                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
8914                          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]);
                         /* GENERATOR SNIP_PDEBC_SYSTEM_2 TOP */  
8915                          ///////////////                          ///////////////
8916                          // process d //                          // process d //
8917                          ///////////////                          ///////////////
# Line 8991  void Brick::assemblePDEBoundarySystem(Pa Line 8956  void Brick::assemblePDEBoundarySystem(Pa
8956                              } else { // constant data                              } else { // constant data
8957                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8958                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8959                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double wd0 = 4*d_p[INDEX2(k, m, numEq)]*w7;
8960                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=16*d_0*w7;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=4*wd0;
8961                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=2*wd0;
8962                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=2*wd0;
8963                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=4*d_0*w7;                                          EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=  wd0;
8964                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=2*wd0;
8965                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=16*d_0*w7;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=4*wd0;
8966                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=4*d_0*w7;                                          EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=  wd0;
8967                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=2*wd0;
8968                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=2*wd0;
8969                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=4*d_0*w7;                                          EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=  wd0;
8970                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=16*d_0*w7;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=4*wd0;
8971                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=2*wd0;
8972                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=4*d_0*w7;                                          EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=  wd0;
8973                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=2*wd0;
8974                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=2*wd0;
8975                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=16*d_0*w7;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=4*wd0;
8976                                      }                                      }
8977                                  }                                  }
8978                              }                              }
# Line 9039  void Brick::assemblePDEBoundarySystem(Pa Line 9004  void Brick::assemblePDEBoundarySystem(Pa
9004                                  }                                  }
9005                              }                              }
9006                          }                          }
                         /* GENERATOR SNIP_PDEBC_SYSTEM_2 BOTTOM */  
9007                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+k0;                          const index_t firstNode=m_NN[0]*m_NN[1]*k2+k0;
9008                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9009                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 9056  void Brick::assemblePDEBoundarySystem(Pa Line 9020  void Brick::assemblePDEBoundarySystem(Pa
9020                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
9021                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
9022                          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]);
                         /* GENERATOR SNIP_PDEBC_SYSTEM_3 TOP */  
9023                          ///////////////                          ///////////////
9024                          // process d //                          // process d //
9025                          ///////////////                          ///////////////
# Line 9101  void Brick::assemblePDEBoundarySystem(Pa Line 9064  void Brick::assemblePDEBoundarySystem(Pa
9064                              } else { // constant data                              } else { // constant data
9065                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9066                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9067                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double wd0 = 4*d_p[INDEX2(k, m, numEq)]*w7;
9068                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=16*d_0*w7;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=4*wd0;
9069                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=2*wd0;
9070                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=2*wd0;
9071                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=4*d_0*w7;                                          EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=  wd0;
9072                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=2*wd0;
9073                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=16*d_0*w7;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=4*wd0;
9074                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=4*d_0*w7;                                          EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=  wd0;
9075                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=2*wd0;
9076                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=2*wd0;
9077                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=4*d_0*w7;                                          EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=  wd0;
9078                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=16*d_0*w7;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=4*wd0;
9079                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=2*wd0;
9080                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=4*d_0*w7;                                          EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=  wd0;
9081                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=2*wd0;
9082                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=8*d_0*w7;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=2*wd0;
9083                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=16*d_0*w7;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=4*wd0;
9084                                      }                                      }
9085                                  }                                  }
9086                              }                              }
# Line 9149  void Brick::assemblePDEBoundarySystem(Pa Line 9112  void Brick::assemblePDEBoundarySystem(Pa
9112                                  }                                  }
9113                              }                              }
9114                          }                          }
                         /* GENERATOR SNIP_PDEBC_SYSTEM_3 BOTTOM */  
9115                          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;
9116                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9117                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 9166  void Brick::assemblePDEBoundarySystem(Pa Line 9128  void Brick::assemblePDEBoundarySystem(Pa
9128                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
9129                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
9130                          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]);
                         /* GENERATOR SNIP_PDEBC_SYSTEM_4 TOP */  
9131                          ///////////////                          ///////////////
9132                          // process d //                          // process d //
9133                          ///////////////                          ///////////////
# Line 9211  void Brick::assemblePDEBoundarySystem(Pa Line 9172  void Brick::assemblePDEBoundarySystem(Pa
9172                              } else { // constant data                              } else { // constant data
9173                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9174                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9175                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double wd0 = 4*d_p[INDEX2(k, m, numEq)]*w12;
9176                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=16*d_0*w12;                                          EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=4*wd0;
9177                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=2*wd0;
9178                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=2*wd0;
9179                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=4*d_0*w12;                                          EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=  wd0;
9180                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=2*wd0;
9181                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=16*d_0*w12;                                          EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=4*wd0;
9182                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=4*d_0*w12;                                          EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=  wd0;
9183                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=2*wd0;
9184                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=2*wd0;
9185                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=4*d_0*w12;                                          EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=  wd0;
9186                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=16*d_0*w12;                                          EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=4*wd0;
9187                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=2*wd0;
9188                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=4*d_0*w12;                                          EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=  wd0;
9189                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=2*wd0;
9190                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=2*wd0;
9191                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=16*d_0*w12;                                          EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=4*wd0;
9192                                      }                                      }
9193                                  }                                  }
9194                              }                              }
# Line 9259  void Brick::assemblePDEBoundarySystem(Pa Line 9220  void Brick::assemblePDEBoundarySystem(Pa
9220                                  }                                  }
9221                              }                              }
9222                          }                          }
                         /* GENERATOR SNIP_PDEBC_SYSTEM_4 BOTTOM */  
9223                          const index_t firstNode=m_NN[0]*k1+k0;                          const index_t firstNode=m_NN[0]*k1+k0;
9224                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9225                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);
# Line 9276  void Brick::assemblePDEBoundarySystem(Pa Line 9236  void Brick::assemblePDEBoundarySystem(Pa
9236                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
9237                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
9238                          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]);
                         /* GENERATOR SNIP_PDEBC_SYSTEM_5 TOP */  
9239                          ///////////////                          ///////////////
9240                          // process d //                          // process d //
9241                          ///////////////                          ///////////////
# Line 9321  void Brick::assemblePDEBoundarySystem(Pa Line 9280  void Brick::assemblePDEBoundarySystem(Pa
9280                              } else { // constant data                              } else { // constant data
9281                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
9282                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
9283                                          const double d_0 = d_p[INDEX2(k, m, numEq)];                                          const double wd0 = 4*d_p[INDEX2(k, m, numEq)]*w12;
9284                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=16*d_0*w12;                                          EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=4*wd0;
9285                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=2*wd0;
9286                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=2*wd0;
9287                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=4*d_0*w12;                                          EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=  wd0;
9288                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=2*wd0;
9289                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=16*d_0*w12;                                          EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=4*wd0;
9290                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=4*d_0*w12;                                          EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=  wd0;
9291                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=2*wd0;
9292                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=2*wd0;
9293                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=4*d_0*w12;                                          EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=  wd0;
9294                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=16*d_0*w12;                                          EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=4*wd0;
9295                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=2*wd0;
9296                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=4*d_0*w12;                                          EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=  wd0;
9297                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=2*wd0;
9298                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=8*d_0*w12;                                          EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=2*wd0;
9299                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=16*d_0*w12;                                          EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=4*wd0;
9300                                      }                                      }
9301                                  }                                  }
9302                              }                              }
# Line 9369  void Brick::assemblePDEBoundarySystem(Pa Line 9328  void Brick::assemblePDEBoundarySystem(Pa
9328                                  }                                  }
9329                              }                              }
9330                          }                          }
                         /* GENERATOR SNIP_PDEBC_SYSTEM_5 BOTTOM */  
9331                          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;
9332                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9333                                  add_EM_F, firstNode, numEq, numComp);                                  add_EM_F, firstNode, numEq, numComp);

Legend:
Removed from v.4378  
changed lines
  Added in v.4382

  ViewVC Help
Powered by ViewVC 1.1.26