/[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 3769 by caltinay, Tue Jan 17 04:09:55 2012 UTC revision 3777 by caltinay, Thu Jan 19 06:17:38 2012 UTC
# Line 617  void Brick::setToSize(escript::Data& out Line 617  void Brick::setToSize(escript::Data& out
617  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
618                                              bool reducedColOrder) const                                              bool reducedColOrder) const
619  {  {
620        /* FIXME: reduced
621      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
622          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
623        */
624      return m_pattern;      return m_pattern;
625  }  }
626    
# Line 1913  void Brick::createPattern() Line 1914  void Brick::createPattern()
1914      Paso_Pattern_free(rowPattern);      Paso_Pattern_free(rowPattern);
1915  }  }
1916    
1917    //private
1918    void Brick::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
1919             const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
1920             bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
1921    {
1922        IndexVector rowIndex;
1923        rowIndex.push_back(m_dofMap[firstNode]);
1924        rowIndex.push_back(m_dofMap[firstNode+1]);
1925        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
1926        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
1927        rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1]);
1928        rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1+1]);
1929        rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);
1930        rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);
1931        if (addF) {
1932            double *F_p=F.getSampleDataRW(0);
1933            for (index_t i=0; i<rowIndex.size(); i++) {
1934                if (rowIndex[i]<getNumDOF()) {
1935                    for (index_t eq=0; eq<nEq; eq++) {
1936                        F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
1937                    }
1938                }
1939            }
1940        }
1941        if (addS) {
1942            addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
1943        }
1944    }
1945    
1946  //protected  //protected
1947  void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in,  void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in,
1948                                         bool reduced) const                                         bool reduced) const
# Line 5286  void Brick::assemblePDESingle(Paso_Syste Line 5316  void Brick::assemblePDESingle(Paso_Syste
5316    
5317                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)
5318                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
5319                          IndexVector rowIndex;                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
5320                          rowIndex.push_back(m_dofMap[firstNode]);                                  add_EM_F, firstNode);
                         rowIndex.push_back(m_dofMap[firstNode+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);  
                         if (add_EM_F) {  
                             double *F_p=rhs.getSampleDataRW(0);  
                             for (index_t i=0; i<rowIndex.size(); i++) {  
                                 if (rowIndex[i]<getNumDOF()) {  
                                     F_p[rowIndex[i]]+=EM_F[i];  
                                 }  
                             }  
                         }  
                         if (add_EM_S) {  
                             addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);  
                         }  
5321                      } // end k0 loop                      } // end k0 loop
5322                  } // end k1 loop                  } // end k1 loop
5323              } // end k2 loop              } // end k2 loop
# Line 5745  void Brick::assemblePDESingleReduced(Pas Line 5757  void Brick::assemblePDESingleReduced(Pas
5757    
5758                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)
5759                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
5760                          IndexVector rowIndex;                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
5761                          rowIndex.push_back(m_dofMap[firstNode]);                                  add_EM_F, firstNode);
                         rowIndex.push_back(m_dofMap[firstNode+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);  
                         if (add_EM_F) {  
                             double *F_p=rhs.getSampleDataRW(0);  
                             for (index_t i=0; i<rowIndex.size(); i++) {  
                                 if (rowIndex[i]<getNumDOF()) {  
                                     F_p[rowIndex[i]]+=EM_F[i];  
                                 }  
                             }  
                         }  
                         if (add_EM_S) {  
                             addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);  
                         }  
5762                      } // end k0 loop                      } // end k0 loop
5763                  } // end k1 loop                  } // end k1 loop
5764              } // end k2 loop              } // end k2 loop
# Line 8911  void Brick::assemblePDESystem(Paso_Syste Line 8905  void Brick::assemblePDESystem(Paso_Syste
8905    
8906                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)
8907                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
8908                          IndexVector rowIndex;                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
8909                          rowIndex.push_back(m_dofMap[firstNode]);                                  add_EM_F, firstNode, numEq, numComp);
                         rowIndex.push_back(m_dofMap[firstNode+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);  
                         if (add_EM_F) {  
                             double *F_p=rhs.getSampleDataRW(0);  
                             for (index_t i=0; i<rowIndex.size(); i++) {  
                                 if (rowIndex[i]<getNumDOF()) {  
                                     for (index_t eq=0; eq<numEq; eq++) {  
                                         F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];  
                                     }  
                                 }  
                             }  
                         }  
                         if (add_EM_S) {  
                             addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);  
                         }  
8910                      } // end k0 loop                      } // end k0 loop
8911                  } // end k1 loop                  } // end k1 loop
8912              } // end k2 loop              } // end k2 loop
# Line 9401  void Brick::assemblePDESystemReduced(Pas Line 9375  void Brick::assemblePDESystemReduced(Pas
9375    
9376                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)
9377                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
9378                          IndexVector rowIndex;                          addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9379                          rowIndex.push_back(m_dofMap[firstNode]);                                  add_EM_F, firstNode, numEq, numComp);
                         rowIndex.push_back(m_dofMap[firstNode+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1+1]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);  
                         rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);  
                         if (add_EM_F) {  
                             double *F_p=rhs.getSampleDataRW(0);  
                             for (index_t i=0; i<rowIndex.size(); i++) {  
                                 if (rowIndex[i]<getNumDOF()) {  
                                     for (index_t eq=0; eq<numEq; eq++) {  
                                         F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];  
                                     }  
                                 }  
                             }  
                         }  
                         if (add_EM_S) {  
                             addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);  
                         }  
9380                      } // end k0 loop                      } // end k0 loop
9381                  } // end k1 loop                  } // end k1 loop
9382              } // end k2 loop              } // end k2 loop
# Line 9434  void Brick::assemblePDESystemReduced(Pas Line 9388  void Brick::assemblePDESystemReduced(Pas
9388  void Brick::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,  void Brick::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
9389        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
9390  {  {
9391        const double h0 = m_l0/m_gNE0;
9392        const double h1 = m_l1/m_gNE1;
9393        const double h2 = m_l2/m_gNE2;
9394        /* GENERATOR SNIP_PDEBC_SINGLE_PRE TOP */
9395        const double w0 = 0.0018607582807716854616*h1*h2;
9396        const double w1 = 0.025917019497006092316*h1*h2;
9397        const double w2 = 0.0069444444444444444444*h1*h2;
9398        const double w3 = 0.00049858867864229740201*h1*h2;
9399        const double w4 = 0.09672363354357992482*h1*h2;
9400        const double w5 = 0.055555555555555555556*h1*h2;
9401        const double w6 = 0.11111111111111111111*h1*h2;
9402        const double w7 = 0.027777777777777777778*h1*h2;
9403        const double w8 = 0.1555021169820365539*h1*h2;
9404        const double w9 = 0.041666666666666666667*h1*h2;
9405        const double w10 = 0.01116454968463011277*h1*h2;
9406        const double w11 = 0.25*h1*h2;
9407        const double w12 = 0.0018607582807716854616*h0*h2;
9408        const double w13 = 0.025917019497006092316*h0*h2;
9409        const double w14 = 0.0069444444444444444444*h0*h2;
9410        const double w15 = 0.00049858867864229740201*h0*h2;
9411        const double w16 = 0.09672363354357992482*h0*h2;
9412        const double w17 = 0.055555555555555555556*h0*h2;
9413        const double w18 = 0.11111111111111111111*h0*h2;
9414        const double w19 = 0.027777777777777777778*h0*h2;
9415        const double w20 = 0.1555021169820365539*h0*h2;
9416        const double w21 = 0.041666666666666666667*h0*h2;
9417        const double w22 = 0.01116454968463011277*h0*h2;
9418        const double w23 = 0.25*h0*h2;
9419        const double w24 = 0.0018607582807716854616*h0*h1;
9420        const double w25 = 0.025917019497006092316*h0*h1;
9421        const double w26 = 0.0069444444444444444444*h0*h1;
9422        const double w27 = 0.00049858867864229740201*h0*h1;
9423        const double w28 = 0.09672363354357992482*h0*h1;
9424        const double w29 = 0.055555555555555555556*h0*h1;
9425        const double w30 = 0.027777777777777777778*h0*h1;
9426        const double w31 = 0.11111111111111111111*h0*h1;
9427        const double w32 = 0.1555021169820365539*h0*h1;
9428        const double w33 = 0.041666666666666666667*h0*h1;
9429        const double w34 = 0.01116454968463011277*h0*h1;
9430        const double w35 = 0.25*h0*h1;
9431        /* GENERATOR SNIP_PDEBC_SINGLE_PRE BOTTOM */
9432        const bool add_EM_S=!d.isEmpty();
9433        const bool add_EM_F=!y.isEmpty();
9434        rhs.requireWrite();
9435    #pragma omp parallel
9436        {
9437            if (m_faceOffset[0] > -1) {
9438                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
9439    #pragma omp for nowait
9440                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
9441                        for (index_t k1=0; k1<m_NE1; ++k1) {
9442                            vector<double> EM_S(8*8, 0);
9443                            vector<double> EM_F(8, 0);
9444                            const index_t e = INDEX2(k1,k2,m_NE1);
9445                            /* GENERATOR SNIP_PDEBC_SINGLE_0 TOP */
9446                            ///////////////
9447                            // process d //
9448                            ///////////////
9449                            if (add_EM_S) {
9450                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
9451                                if (d.actsExpanded()) {
9452                                    const double d_0 = d_p[0];
9453                                    const double d_1 = d_p[1];
9454                                    const double d_2 = d_p[2];
9455                                    const double d_3 = d_p[3];
9456                                    const double tmp3_0 = d_1 + d_3;
9457                                    const double tmp6_0 = d_0 + d_1 + d_2 + d_3;
9458                                    const double tmp2_0 = d_0 + d_2;
9459                                    const double tmp0_0 = d_0 + d_1;
9460                                    const double tmp4_0 = d_0 + d_3;
9461                                    const double tmp1_0 = d_2 + d_3;
9462                                    const double tmp5_0 = d_1 + d_2;
9463                                    const double tmp14_1 = d_0*w4;
9464                                    const double tmp17_1 = d_3*w4;
9465                                    const double tmp12_1 = d_1*w4;
9466                                    const double tmp9_1 = d_2*w4;
9467                                    const double tmp4_1 = tmp3_0*w0;
9468                                    const double tmp3_1 = tmp3_0*w1;
9469                                    const double tmp7_1 = tmp0_0*w1;
9470                                    const double tmp8_1 = d_1*w3;
9471                                    const double tmp16_1 = d_0*w3;
9472                                    const double tmp18_1 = tmp6_0*w2;
9473                                    const double tmp1_1 = tmp1_0*w1;
9474                                    const double tmp15_1 = tmp5_0*w2;
9475                                    const double tmp6_1 = tmp1_0*w0;
9476                                    const double tmp2_1 = tmp2_0*w0;
9477                                    const double tmp13_1 = d_3*w3;
9478                                    const double tmp5_1 = tmp2_0*w1;
9479                                    const double tmp10_1 = tmp4_0*w2;
9480                                    const double tmp11_1 = d_2*w3;
9481                                    const double tmp0_1 = tmp0_0*w0;
9482                                    EM_S[INDEX2(6,4,8)]+=tmp0_1 + tmp1_1;
9483                                    EM_S[INDEX2(2,6,8)]+=tmp2_1 + tmp3_1;
9484                                    EM_S[INDEX2(4,0,8)]+=tmp4_1 + tmp5_1;
9485                                    EM_S[INDEX2(2,0,8)]+=tmp6_1 + tmp7_1;
9486                                    EM_S[INDEX2(4,4,8)]+=tmp10_1 + tmp8_1 + tmp9_1;
9487                                    EM_S[INDEX2(2,2,8)]+=tmp10_1 + tmp11_1 + tmp12_1;
9488                                    EM_S[INDEX2(0,0,8)]+=tmp13_1 + tmp14_1 + tmp15_1;
9489                                    EM_S[INDEX2(6,6,8)]+=tmp15_1 + tmp16_1 + tmp17_1;
9490                                    EM_S[INDEX2(0,4,8)]+=tmp4_1 + tmp5_1;
9491                                    EM_S[INDEX2(6,0,8)]+=tmp18_1;
9492                                    EM_S[INDEX2(4,2,8)]+=tmp18_1;
9493                                    EM_S[INDEX2(4,6,8)]+=tmp0_1 + tmp1_1;
9494                                    EM_S[INDEX2(0,2,8)]+=tmp6_1 + tmp7_1;
9495                                    EM_S[INDEX2(0,6,8)]+=tmp18_1;
9496                                    EM_S[INDEX2(6,2,8)]+=tmp2_1 + tmp3_1;
9497                                    EM_S[INDEX2(2,4,8)]+=tmp18_1;
9498                                } else { /* constant data */
9499                                    const double tmp0_1 = d_p[0]*w5;
9500                                    const double tmp2_1 = d_p[0]*w7;
9501                                    const double tmp1_1 = d_p[0]*w6;
9502                                    EM_S[INDEX2(6,4,8)]+=tmp0_1;
9503                                    EM_S[INDEX2(2,6,8)]+=tmp0_1;
9504                                    EM_S[INDEX2(4,0,8)]+=tmp0_1;
9505                                    EM_S[INDEX2(2,0,8)]+=tmp0_1;
9506                                    EM_S[INDEX2(4,4,8)]+=tmp1_1;
9507                                    EM_S[INDEX2(2,2,8)]+=tmp1_1;
9508                                    EM_S[INDEX2(0,0,8)]+=tmp1_1;
9509                                    EM_S[INDEX2(6,6,8)]+=tmp1_1;
9510                                    EM_S[INDEX2(0,4,8)]+=tmp0_1;
9511                                    EM_S[INDEX2(6,0,8)]+=tmp2_1;
9512                                    EM_S[INDEX2(4,2,8)]+=tmp2_1;
9513                                    EM_S[INDEX2(4,6,8)]+=tmp0_1;
9514                                    EM_S[INDEX2(0,2,8)]+=tmp0_1;
9515                                    EM_S[INDEX2(0,6,8)]+=tmp2_1;
9516                                    EM_S[INDEX2(6,2,8)]+=tmp0_1;
9517                                    EM_S[INDEX2(2,4,8)]+=tmp2_1;
9518                                }
9519                            }
9520                            ///////////////
9521                            // process y //
9522                            ///////////////
9523                            if (add_EM_F) {
9524                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
9525                                if (y.actsExpanded()) {
9526                                    const double y_0 = y_p[0];
9527                                    const double y_1 = y_p[1];
9528                                    const double y_2 = y_p[2];
9529                                    const double y_3 = y_p[3];
9530                                    const double tmp0_0 = y_1 + y_2;
9531                                    const double tmp1_0 = y_0 + y_3;
9532                                    const double tmp2_1 = w10*y_3;
9533                                    const double tmp8_1 = w10*y_0;
9534                                    const double tmp5_1 = w10*y_2;
9535                                    const double tmp3_1 = w8*y_1;
9536                                    const double tmp9_1 = w8*y_3;
9537                                    const double tmp0_1 = w8*y_0;
9538                                    const double tmp1_1 = tmp0_0*w9;
9539                                    const double tmp7_1 = w8*y_2;
9540                                    const double tmp4_1 = tmp1_0*w9;
9541                                    const double tmp6_1 = w10*y_1;
9542                                    EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
9543                                    EM_F[2]+=tmp3_1 + tmp4_1 + tmp5_1;
9544                                    EM_F[4]+=tmp4_1 + tmp6_1 + tmp7_1;
9545                                    EM_F[6]+=tmp1_1 + tmp8_1 + tmp9_1;
9546                                } else { /* constant data */
9547                                    const double tmp0_1 = w11*y_p[0];
9548                                    EM_F[0]+=tmp0_1;
9549                                    EM_F[2]+=tmp0_1;
9550                                    EM_F[4]+=tmp0_1;
9551                                    EM_F[6]+=tmp0_1;
9552                                }
9553                            }
9554                            /* GENERATOR SNIP_PDEBC_SINGLE_0 BOTTOM */
9555                            const index_t firstNode=m_N0*m_N1*k2+m_N0*k1;
9556                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9557                                    add_EM_F, firstNode);
9558                        } // k1 loop
9559                    } // k2 loop
9560                } // colouring
9561            } // face 0
9562    
9563            if (m_faceOffset[1] > -1) {
9564                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
9565    #pragma omp for nowait
9566                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
9567                        for (index_t k1=0; k1<m_NE1; ++k1) {
9568                            vector<double> EM_S(8*8, 0);
9569                            vector<double> EM_F(8, 0);
9570                            const index_t e = m_faceOffset[1]+INDEX2(k1,k2,m_NE1);
9571                            /* GENERATOR SNIP_PDEBC_SINGLE_1 TOP */
9572                            ///////////////
9573                            // process d //
9574                            ///////////////
9575                            if (add_EM_S) {
9576                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
9577                                if (d.actsExpanded()) {
9578                                    const double d_0 = d_p[0];
9579                                    const double d_1 = d_p[1];
9580                                    const double d_2 = d_p[2];
9581                                    const double d_3 = d_p[3];
9582                                    const double tmp1_0 = d_1 + d_3;
9583                                    const double tmp6_0 = d_0 + d_1 + d_2 + d_3;
9584                                    const double tmp0_0 = d_0 + d_2;
9585                                    const double tmp3_0 = d_0 + d_1;
9586                                    const double tmp4_0 = d_0 + d_3;
9587                                    const double tmp2_0 = d_2 + d_3;
9588                                    const double tmp5_0 = d_1 + d_2;
9589                                    const double tmp10_1 = d_3*w4;
9590                                    const double tmp13_1 = tmp2_0*w1;
9591                                    const double tmp16_1 = d_0*w4;
9592                                    const double tmp18_1 = d_2*w4;
9593                                    const double tmp12_1 = tmp3_0*w0;
9594                                    const double tmp3_1 = tmp3_0*w1;
9595                                    const double tmp5_1 = tmp0_0*w1;
9596                                    const double tmp17_1 = d_1*w3;
9597                                    const double tmp9_1 = d_0*w3;
9598                                    const double tmp14_1 = tmp6_0*w2;
9599                                    const double tmp1_1 = tmp1_0*w1;
9600                                    const double tmp11_1 = tmp5_0*w2;
9601                                    const double tmp4_1 = tmp1_0*w0;
9602                                    const double tmp2_1 = tmp2_0*w0;
9603                                    const double tmp15_1 = d_3*w3;
9604                                    const double tmp7_1 = d_1*w4;
9605                                    const double tmp8_1 = tmp4_0*w2;
9606                                    const double tmp6_1 = d_2*w3;
9607                                    const double tmp0_1 = tmp0_0*w0;
9608                                    EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1;
9609                                    EM_S[INDEX2(1,3,8)]+=tmp2_1 + tmp3_1;
9610                                    EM_S[INDEX2(5,1,8)]+=tmp4_1 + tmp5_1;
9611                                    EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp1_1;
9612                                    EM_S[INDEX2(3,3,8)]+=tmp6_1 + tmp7_1 + tmp8_1;
9613                                    EM_S[INDEX2(1,5,8)]+=tmp4_1 + tmp5_1;
9614                                    EM_S[INDEX2(7,7,8)]+=tmp10_1 + tmp11_1 + tmp9_1;
9615                                    EM_S[INDEX2(5,7,8)]+=tmp12_1 + tmp13_1;
9616                                    EM_S[INDEX2(5,3,8)]+=tmp14_1;
9617                                    EM_S[INDEX2(1,1,8)]+=tmp11_1 + tmp15_1 + tmp16_1;
9618                                    EM_S[INDEX2(7,1,8)]+=tmp14_1;
9619                                    EM_S[INDEX2(5,5,8)]+=tmp17_1 + tmp18_1 + tmp8_1;
9620                                    EM_S[INDEX2(7,5,8)]+=tmp12_1 + tmp13_1;
9621                                    EM_S[INDEX2(3,5,8)]+=tmp14_1;
9622                                    EM_S[INDEX2(3,1,8)]+=tmp2_1 + tmp3_1;
9623                                    EM_S[INDEX2(1,7,8)]+=tmp14_1;
9624                                } else { /* constant data */
9625                                    const double tmp0_1 = d_p[0]*w5;
9626                                    const double tmp2_1 = d_p[0]*w7;
9627                                    const double tmp1_1 = d_p[0]*w6;
9628                                    EM_S[INDEX2(7,3,8)]+=tmp0_1;
9629                                    EM_S[INDEX2(1,3,8)]+=tmp0_1;
9630                                    EM_S[INDEX2(5,1,8)]+=tmp0_1;
9631                                    EM_S[INDEX2(3,7,8)]+=tmp0_1;
9632                                    EM_S[INDEX2(3,3,8)]+=tmp1_1;
9633                                    EM_S[INDEX2(1,5,8)]+=tmp0_1;
9634                                    EM_S[INDEX2(7,7,8)]+=tmp1_1;
9635                                    EM_S[INDEX2(5,7,8)]+=tmp0_1;
9636                                    EM_S[INDEX2(5,3,8)]+=tmp2_1;
9637                                    EM_S[INDEX2(1,1,8)]+=tmp1_1;
9638                                    EM_S[INDEX2(7,1,8)]+=tmp2_1;
9639                                    EM_S[INDEX2(5,5,8)]+=tmp1_1;
9640                                    EM_S[INDEX2(7,5,8)]+=tmp0_1;
9641                                    EM_S[INDEX2(3,5,8)]+=tmp2_1;
9642                                    EM_S[INDEX2(3,1,8)]+=tmp0_1;
9643                                    EM_S[INDEX2(1,7,8)]+=tmp2_1;
9644                                }
9645                            }
9646                            ///////////////
9647                            // process y //
9648                            ///////////////
9649                            if (add_EM_F) {
9650                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
9651                                if (y.actsExpanded()) {
9652                                    const double y_0 = y_p[0];
9653                                    const double y_1 = y_p[1];
9654                                    const double y_2 = y_p[2];
9655                                    const double y_3 = y_p[3];
9656                                    const double tmp0_0 = y_1 + y_2;
9657                                    const double tmp1_0 = y_0 + y_3;
9658                                    const double tmp2_1 = w10*y_3;
9659                                    const double tmp8_1 = w10*y_0;
9660                                    const double tmp5_1 = w10*y_2;
9661                                    const double tmp3_1 = w8*y_1;
9662                                    const double tmp9_1 = w8*y_3;
9663                                    const double tmp0_1 = w8*y_0;
9664                                    const double tmp1_1 = tmp0_0*w9;
9665                                    const double tmp7_1 = w8*y_2;
9666                                    const double tmp4_1 = tmp1_0*w9;
9667                                    const double tmp6_1 = w10*y_1;
9668                                    EM_F[1]+=tmp0_1 + tmp1_1 + tmp2_1;
9669                                    EM_F[3]+=tmp3_1 + tmp4_1 + tmp5_1;
9670                                    EM_F[5]+=tmp4_1 + tmp6_1 + tmp7_1;
9671                                    EM_F[7]+=tmp1_1 + tmp8_1 + tmp9_1;
9672                                } else { /* constant data */
9673                                    const double tmp0_1 = w11*y_p[0];
9674                                    EM_F[1]+=tmp0_1;
9675                                    EM_F[3]+=tmp0_1;
9676                                    EM_F[5]+=tmp0_1;
9677                                    EM_F[7]+=tmp0_1;
9678                                }
9679                            }
9680                            /* GENERATOR SNIP_PDEBC_SINGLE_1 BOTTOM */
9681                            const index_t firstNode=m_N0*m_N1*k2+m_N0*(k1+1)-2;
9682                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9683                                    add_EM_F, firstNode);
9684                        } // k1 loop
9685                    } // k2 loop
9686                } // colouring
9687            } // face 1
9688    
9689            if (m_faceOffset[2] > -1) {
9690                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
9691    #pragma omp for nowait
9692                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
9693                        for (index_t k0=0; k0<m_NE0; ++k0) {
9694                            vector<double> EM_S(8*8, 0);
9695                            vector<double> EM_F(8, 0);
9696                            const index_t e = m_faceOffset[2]+INDEX2(k0,k2,m_NE0);
9697                            /* GENERATOR SNIP_PDEBC_SINGLE_2 TOP */
9698                            ///////////////
9699                            // process d //
9700                            ///////////////
9701                            if (add_EM_S) {
9702                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
9703                                if (d.actsExpanded()) {
9704                                    const double d_0 = d_p[0];
9705                                    const double d_1 = d_p[1];
9706                                    const double d_2 = d_p[2];
9707                                    const double d_3 = d_p[3];
9708                                    const double tmp2_0 = d_1 + d_3;
9709                                    const double tmp5_0 = d_0 + d_1 + d_2 + d_3;
9710                                    const double tmp3_0 = d_0 + d_2;
9711                                    const double tmp1_0 = d_0 + d_1;
9712                                    const double tmp4_0 = d_0 + d_3;
9713                                    const double tmp0_0 = d_2 + d_3;
9714                                    const double tmp6_0 = d_1 + d_2;
9715                                    const double tmp2_1 = tmp2_0*w13;
9716                                    const double tmp14_1 = d_3*w15;
9717                                    const double tmp0_1 = tmp0_0*w13;
9718                                    const double tmp3_1 = tmp3_0*w12;
9719                                    const double tmp17_1 = tmp1_0*w13;
9720                                    const double tmp18_1 = tmp0_0*w12;
9721                                    const double tmp8_1 = d_1*w15;
9722                                    const double tmp16_1 = d_0*w15;
9723                                    const double tmp11_1 = d_2*w15;
9724                                    const double tmp5_1 = tmp2_0*w12;
9725                                    const double tmp15_1 = d_3*w16;
9726                                    const double tmp13_1 = tmp6_0*w14;
9727                                    const double tmp1_1 = tmp1_0*w12;
9728                                    const double tmp7_1 = tmp4_0*w14;
9729                                    const double tmp12_1 = d_0*w16;
9730                                    const double tmp10_1 = d_1*w16;
9731                                    const double tmp6_1 = d_2*w16;
9732                                    const double tmp9_1 = tmp5_0*w14;
9733                                    const double tmp4_1 = tmp3_0*w13;
9734                                    EM_S[INDEX2(5,4,8)]+=tmp0_1 + tmp1_1;
9735                                    EM_S[INDEX2(5,1,8)]+=tmp2_1 + tmp3_1;
9736                                    EM_S[INDEX2(4,0,8)]+=tmp4_1 + tmp5_1;
9737                                    EM_S[INDEX2(4,4,8)]+=tmp6_1 + tmp7_1 + tmp8_1;
9738                                    EM_S[INDEX2(1,5,8)]+=tmp2_1 + tmp3_1;
9739                                    EM_S[INDEX2(4,1,8)]+=tmp9_1;
9740                                    EM_S[INDEX2(1,1,8)]+=tmp10_1 + tmp11_1 + tmp7_1;
9741                                    EM_S[INDEX2(0,0,8)]+=tmp12_1 + tmp13_1 + tmp14_1;
9742                                    EM_S[INDEX2(5,0,8)]+=tmp9_1;
9743                                    EM_S[INDEX2(4,5,8)]+=tmp0_1 + tmp1_1;
9744                                    EM_S[INDEX2(0,4,8)]+=tmp4_1 + tmp5_1;
9745                                    EM_S[INDEX2(5,5,8)]+=tmp13_1 + tmp15_1 + tmp16_1;
9746                                    EM_S[INDEX2(1,4,8)]+=tmp9_1;
9747                                    EM_S[INDEX2(1,0,8)]+=tmp17_1 + tmp18_1;
9748                                    EM_S[INDEX2(0,1,8)]+=tmp17_1 + tmp18_1;
9749                                    EM_S[INDEX2(0,5,8)]+=tmp9_1;
9750                                } else { /* constant data */
9751                                    const double tmp0_1 = d_p[0]*w17;
9752                                    const double tmp2_1 = d_p[0]*w19;
9753                                    const double tmp1_1 = d_p[0]*w18;
9754                                    EM_S[INDEX2(5,4,8)]+=tmp0_1;
9755                                    EM_S[INDEX2(5,1,8)]+=tmp0_1;
9756                                    EM_S[INDEX2(4,0,8)]+=tmp0_1;
9757                                    EM_S[INDEX2(4,4,8)]+=tmp1_1;
9758                                    EM_S[INDEX2(1,5,8)]+=tmp0_1;
9759                                    EM_S[INDEX2(4,1,8)]+=tmp2_1;
9760                                    EM_S[INDEX2(1,1,8)]+=tmp1_1;
9761                                    EM_S[INDEX2(0,0,8)]+=tmp1_1;
9762                                    EM_S[INDEX2(5,0,8)]+=tmp2_1;
9763                                    EM_S[INDEX2(4,5,8)]+=tmp0_1;
9764                                    EM_S[INDEX2(0,4,8)]+=tmp0_1;
9765                                    EM_S[INDEX2(5,5,8)]+=tmp1_1;
9766                                    EM_S[INDEX2(1,4,8)]+=tmp2_1;
9767                                    EM_S[INDEX2(1,0,8)]+=tmp0_1;
9768                                    EM_S[INDEX2(0,1,8)]+=tmp0_1;
9769                                    EM_S[INDEX2(0,5,8)]+=tmp2_1;
9770                                }
9771                            }
9772                            ///////////////
9773                            // process y //
9774                            ///////////////
9775                            if (add_EM_F) {
9776                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
9777                                if (y.actsExpanded()) {
9778                                    const double y_0 = y_p[0];
9779                                    const double y_1 = y_p[1];
9780                                    const double y_2 = y_p[2];
9781                                    const double y_3 = y_p[3];
9782                                    const double tmp0_0 = y_1 + y_2;
9783                                    const double tmp1_0 = y_0 + y_3;
9784                                    const double tmp0_1 = w22*y_3;
9785                                    const double tmp6_1 = w22*y_1;
9786                                    const double tmp3_1 = w22*y_2;
9787                                    const double tmp5_1 = w20*y_1;
9788                                    const double tmp9_1 = w20*y_3;
9789                                    const double tmp4_1 = tmp1_0*w21;
9790                                    const double tmp8_1 = w22*y_0;
9791                                    const double tmp2_1 = w20*y_0;
9792                                    const double tmp7_1 = w20*y_2;
9793                                    const double tmp1_1 = tmp0_0*w21;
9794                                    EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
9795                                    EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
9796                                    EM_F[4]+=tmp4_1 + tmp6_1 + tmp7_1;
9797                                    EM_F[5]+=tmp1_1 + tmp8_1 + tmp9_1;
9798                                } else { /* constant data */
9799                                    const double tmp0_1 = w23*y_p[0];
9800                                    EM_F[0]+=tmp0_1;
9801                                    EM_F[1]+=tmp0_1;
9802                                    EM_F[4]+=tmp0_1;
9803                                    EM_F[5]+=tmp0_1;
9804                                }
9805                            }
9806                            /* GENERATOR SNIP_PDEBC_SINGLE_2 BOTTOM */
9807                            const index_t firstNode=m_N0*m_N1*k2+k0;
9808                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9809                                    add_EM_F, firstNode);
9810                        } // k0 loop
9811                    } // k2 loop
9812                } // colouring
9813            } // face 2
9814    
9815            if (m_faceOffset[3] > -1) {
9816                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
9817    #pragma omp for nowait
9818                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
9819                        for (index_t k0=0; k0<m_NE0; ++k0) {
9820                            vector<double> EM_S(8*8, 0);
9821                            vector<double> EM_F(8, 0);
9822                            const index_t e = m_faceOffset[3]+INDEX2(k0,k2,m_NE0);
9823                            /* GENERATOR SNIP_PDEBC_SINGLE_3 TOP */
9824                            ///////////////
9825                            // process d //
9826                            ///////////////
9827                            if (add_EM_S) {
9828                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
9829                                if (d.actsExpanded()) {
9830                                    const double d_0 = d_p[0];
9831                                    const double d_1 = d_p[1];
9832                                    const double d_2 = d_p[2];
9833                                    const double d_3 = d_p[3];
9834                                    const double tmp0_0 = d_1 + d_3;
9835                                    const double tmp2_0 = d_0 + d_1 + d_2 + d_3;
9836                                    const double tmp1_0 = d_0 + d_2;
9837                                    const double tmp4_0 = d_0 + d_1;
9838                                    const double tmp5_0 = d_0 + d_3;
9839                                    const double tmp3_0 = d_2 + d_3;
9840                                    const double tmp6_0 = d_1 + d_2;
9841                                    const double tmp15_1 = tmp4_0*w13;
9842                                    const double tmp10_1 = d_0*w16;
9843                                    const double tmp6_1 = tmp4_0*w12;
9844                                    const double tmp16_1 = tmp3_0*w12;
9845                                    const double tmp0_1 = tmp0_0*w13;
9846                                    const double tmp2_1 = tmp1_0*w13;
9847                                    const double tmp18_1 = d_1*w15;
9848                                    const double tmp14_1 = d_0*w15;
9849                                    const double tmp9_1 = d_2*w15;
9850                                    const double tmp4_1 = tmp2_0*w14;
9851                                    const double tmp13_1 = d_3*w16;
9852                                    const double tmp11_1 = tmp6_0*w14;
9853                                    const double tmp1_1 = tmp1_0*w12;
9854                                    const double tmp12_1 = d_3*w15;
9855                                    const double tmp3_1 = tmp0_0*w12;
9856                                    const double tmp7_1 = d_1*w16;
9857                                    const double tmp17_1 = d_2*w16;
9858                                    const double tmp8_1 = tmp5_0*w14;
9859                                    const double tmp5_1 = tmp3_0*w13;
9860                                    EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1;
9861                                    EM_S[INDEX2(2,6,8)]+=tmp2_1 + tmp3_1;
9862                                    EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp1_1;
9863                                    EM_S[INDEX2(7,2,8)]+=tmp4_1;
9864                                    EM_S[INDEX2(6,7,8)]+=tmp5_1 + tmp6_1;
9865                                    EM_S[INDEX2(3,3,8)]+=tmp7_1 + tmp8_1 + tmp9_1;
9866                                    EM_S[INDEX2(7,6,8)]+=tmp5_1 + tmp6_1;
9867                                    EM_S[INDEX2(6,3,8)]+=tmp4_1;
9868                                    EM_S[INDEX2(3,6,8)]+=tmp4_1;
9869                                    EM_S[INDEX2(2,2,8)]+=tmp10_1 + tmp11_1 + tmp12_1;
9870                                    EM_S[INDEX2(7,7,8)]+=tmp11_1 + tmp13_1 + tmp14_1;
9871                                    EM_S[INDEX2(2,7,8)]+=tmp4_1;
9872                                    EM_S[INDEX2(3,2,8)]+=tmp15_1 + tmp16_1;
9873                                    EM_S[INDEX2(6,6,8)]+=tmp17_1 + tmp18_1 + tmp8_1;
9874                                    EM_S[INDEX2(2,3,8)]+=tmp15_1 + tmp16_1;
9875                                    EM_S[INDEX2(6,2,8)]+=tmp2_1 + tmp3_1;
9876                                } else { /* constant data */
9877                                    const double tmp0_1 = d_p[0]*w17;
9878                                    const double tmp1_1 = d_p[0]*w19;
9879                                    const double tmp2_1 = d_p[0]*w18;
9880                                    EM_S[INDEX2(7,3,8)]+=tmp0_1;
9881                                    EM_S[INDEX2(2,6,8)]+=tmp0_1;
9882                                    EM_S[INDEX2(3,7,8)]+=tmp0_1;
9883                                    EM_S[INDEX2(7,2,8)]+=tmp1_1;
9884                                    EM_S[INDEX2(6,7,8)]+=tmp0_1;
9885                                    EM_S[INDEX2(3,3,8)]+=tmp2_1;
9886                                    EM_S[INDEX2(7,6,8)]+=tmp0_1;
9887                                    EM_S[INDEX2(6,3,8)]+=tmp1_1;
9888                                    EM_S[INDEX2(3,6,8)]+=tmp1_1;
9889                                    EM_S[INDEX2(2,2,8)]+=tmp2_1;
9890                                    EM_S[INDEX2(7,7,8)]+=tmp2_1;
9891                                    EM_S[INDEX2(2,7,8)]+=tmp1_1;
9892                                    EM_S[INDEX2(3,2,8)]+=tmp0_1;
9893                                    EM_S[INDEX2(6,6,8)]+=tmp2_1;
9894                                    EM_S[INDEX2(2,3,8)]+=tmp0_1;
9895                                    EM_S[INDEX2(6,2,8)]+=tmp0_1;
9896                                }
9897                            }
9898                            ///////////////
9899                            // process y //
9900                            ///////////////
9901                            if (add_EM_F) {
9902                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
9903                                if (y.actsExpanded()) {
9904                                    const double y_0 = y_p[0];
9905                                    const double y_1 = y_p[1];
9906                                    const double y_2 = y_p[2];
9907                                    const double y_3 = y_p[3];
9908                                    const double tmp0_0 = y_1 + y_2;
9909                                    const double tmp1_0 = y_0 + y_3;
9910                                    const double tmp0_1 = w22*y_3;
9911                                    const double tmp6_1 = w22*y_1;
9912                                    const double tmp3_1 = w22*y_2;
9913                                    const double tmp5_1 = w20*y_1;
9914                                    const double tmp9_1 = w20*y_3;
9915                                    const double tmp4_1 = tmp1_0*w21;
9916                                    const double tmp8_1 = w22*y_0;
9917                                    const double tmp2_1 = w20*y_0;
9918                                    const double tmp7_1 = w20*y_2;
9919                                    const double tmp1_1 = tmp0_0*w21;
9920                                    EM_F[2]+=tmp0_1 + tmp1_1 + tmp2_1;
9921                                    EM_F[3]+=tmp3_1 + tmp4_1 + tmp5_1;
9922                                    EM_F[6]+=tmp4_1 + tmp6_1 + tmp7_1;
9923                                    EM_F[7]+=tmp1_1 + tmp8_1 + tmp9_1;
9924                                } else { /* constant data */
9925                                    const double tmp0_1 = w23*y_p[0];
9926                                    EM_F[2]+=tmp0_1;
9927                                    EM_F[3]+=tmp0_1;
9928                                    EM_F[6]+=tmp0_1;
9929                                    EM_F[7]+=tmp0_1;
9930                                }
9931                            }
9932                            /* GENERATOR SNIP_PDEBC_SINGLE_3 BOTTOM */
9933                            const index_t firstNode=m_N0*m_N1*k2+m_N0*(m_N1-2)+k0;
9934                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
9935                                    add_EM_F, firstNode);
9936                        } // k0 loop
9937                    } // k2 loop
9938                } // colouring
9939            } // face 3
9940    
9941            if (m_faceOffset[4] > -1) {
9942                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
9943    #pragma omp for nowait
9944                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
9945                        for (index_t k0=0; k0<m_NE0; ++k0) {
9946                            vector<double> EM_S(8*8, 0);
9947                            vector<double> EM_F(8, 0);
9948                            const index_t e = m_faceOffset[4]+INDEX2(k0,k1,m_NE0);
9949                            /* GENERATOR SNIP_PDEBC_SINGLE_4 TOP */
9950                            ///////////////
9951                            // process d //
9952                            ///////////////
9953                            if (add_EM_S) {
9954                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
9955                                if (d.actsExpanded()) {
9956                                    const double d_0 = d_p[0];
9957                                    const double d_1 = d_p[1];
9958                                    const double d_2 = d_p[2];
9959                                    const double d_3 = d_p[3];
9960                                    const double tmp0_0 = d_1 + d_3;
9961                                    const double tmp2_0 = d_0 + d_1 + d_2 + d_3;
9962                                    const double tmp1_0 = d_0 + d_2;
9963                                    const double tmp6_0 = d_0 + d_1;
9964                                    const double tmp4_0 = d_0 + d_3;
9965                                    const double tmp5_0 = d_2 + d_3;
9966                                    const double tmp3_0 = d_1 + d_2;
9967                                    const double tmp18_1 = tmp5_0*w24;
9968                                    const double tmp6_1 = tmp1_0*w25;
9969                                    const double tmp4_1 = d_0*w27;
9970                                    const double tmp12_1 = d_2*w27;
9971                                    const double tmp0_1 = tmp0_0*w25;
9972                                    const double tmp5_1 = tmp3_0*w26;
9973                                    const double tmp2_1 = tmp2_0*w26;
9974                                    const double tmp17_1 = tmp6_0*w25;
9975                                    const double tmp14_1 = tmp6_0*w24;
9976                                    const double tmp11_1 = d_1*w28;
9977                                    const double tmp9_1 = d_1*w27;
9978                                    const double tmp16_1 = d_3*w27;
9979                                    const double tmp8_1 = d_2*w28;
9980                                    const double tmp7_1 = tmp0_0*w24;
9981                                    const double tmp15_1 = d_0*w28;
9982                                    const double tmp13_1 = tmp5_0*w25;
9983                                    const double tmp3_1 = d_3*w28;
9984                                    const double tmp10_1 = tmp4_0*w26;
9985                                    const double tmp1_1 = tmp1_0*w24;
9986                                    EM_S[INDEX2(1,3,8)]+=tmp0_1 + tmp1_1;
9987                                    EM_S[INDEX2(3,0,8)]+=tmp2_1;
9988                                    EM_S[INDEX2(0,3,8)]+=tmp2_1;
9989                                    EM_S[INDEX2(1,2,8)]+=tmp2_1;
9990                                    EM_S[INDEX2(3,3,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9991                                    EM_S[INDEX2(2,0,8)]+=tmp6_1 + tmp7_1;
9992                                    EM_S[INDEX2(2,2,8)]+=tmp10_1 + tmp8_1 + tmp9_1;
9993                                    EM_S[INDEX2(1,1,8)]+=tmp10_1 + tmp11_1 + tmp12_1;
9994                                    EM_S[INDEX2(3,2,8)]+=tmp13_1 + tmp14_1;
9995                                    EM_S[INDEX2(0,0,8)]+=tmp15_1 + tmp16_1 + tmp5_1;
9996                                    EM_S[INDEX2(2,3,8)]+=tmp13_1 + tmp14_1;
9997                                    EM_S[INDEX2(2,1,8)]+=tmp2_1;
9998                                    EM_S[INDEX2(1,0,8)]+=tmp17_1 + tmp18_1;
9999                                    EM_S[INDEX2(0,1,8)]+=tmp17_1 + tmp18_1;
10000                                    EM_S[INDEX2(3,1,8)]+=tmp0_1 + tmp1_1;
10001                                    EM_S[INDEX2(0,2,8)]+=tmp6_1 + tmp7_1;
10002                                } else { /* constant data */
10003                                    const double tmp2_1 = d_p[0]*w31;
10004                                    const double tmp1_1 = d_p[0]*w30;
10005                                    const double tmp0_1 = d_p[0]*w29;
10006                                    EM_S[INDEX2(1,3,8)]+=tmp0_1;
10007                                    EM_S[INDEX2(3,0,8)]+=tmp1_1;
10008                                    EM_S[INDEX2(0,3,8)]+=tmp1_1;
10009                                    EM_S[INDEX2(1,2,8)]+=tmp1_1;
10010                                    EM_S[INDEX2(3,3,8)]+=tmp2_1;
10011                                    EM_S[INDEX2(2,0,8)]+=tmp0_1;
10012                                    EM_S[INDEX2(2,2,8)]+=tmp2_1;
10013                                    EM_S[INDEX2(1,1,8)]+=tmp2_1;
10014                                    EM_S[INDEX2(3,2,8)]+=tmp0_1;
10015                                    EM_S[INDEX2(0,0,8)]+=tmp2_1;
10016                                    EM_S[INDEX2(2,3,8)]+=tmp0_1;
10017                                    EM_S[INDEX2(2,1,8)]+=tmp1_1;
10018                                    EM_S[INDEX2(1,0,8)]+=tmp0_1;
10019                                    EM_S[INDEX2(0,1,8)]+=tmp0_1;
10020                                    EM_S[INDEX2(3,1,8)]+=tmp0_1;
10021                                    EM_S[INDEX2(0,2,8)]+=tmp0_1;
10022                                }
10023                            }
10024                            ///////////////
10025                            // process y //
10026                            ///////////////
10027                            if (add_EM_F) {
10028                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10029                                if (y.actsExpanded()) {
10030                                    const double y_0 = y_p[0];
10031                                    const double y_1 = y_p[1];
10032                                    const double y_2 = y_p[2];
10033                                    const double y_3 = y_p[3];
10034                                    const double tmp0_0 = y_1 + y_2;
10035                                    const double tmp1_0 = y_0 + y_3;
10036                                    const double tmp7_1 = w34*y_1;
10037                                    const double tmp3_1 = w32*y_1;
10038                                    const double tmp8_1 = w32*y_3;
10039                                    const double tmp0_1 = w32*y_0;
10040                                    const double tmp2_1 = w34*y_3;
10041                                    const double tmp9_1 = w34*y_0;
10042                                    const double tmp6_1 = w32*y_2;
10043                                    const double tmp5_1 = w34*y_2;
10044                                    const double tmp1_1 = tmp0_0*w33;
10045                                    const double tmp4_1 = tmp1_0*w33;
10046                                    EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
10047                                    EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
10048                                    EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
10049                                    EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
10050                                } else { /* constant data */
10051                                    const double tmp0_1 = w35*y_p[0];
10052                                    EM_F[0]+=tmp0_1;
10053                                    EM_F[1]+=tmp0_1;
10054                                    EM_F[2]+=tmp0_1;
10055                                    EM_F[3]+=tmp0_1;
10056                                }
10057                            }
10058                            /* GENERATOR SNIP_PDEBC_SINGLE_4 BOTTOM */
10059                            const index_t firstNode=m_N0*k1+k0;
10060                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10061                                    add_EM_F, firstNode);
10062                        } // k0 loop
10063                    } // k1 loop
10064                } // colouring
10065            } // face 4
10066    
10067            if (m_faceOffset[5] > -1) {
10068                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
10069    #pragma omp for nowait
10070                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
10071                        for (index_t k0=0; k0<m_NE0; ++k0) {
10072                            vector<double> EM_S(8*8, 0);
10073                            vector<double> EM_F(8, 0);
10074                            const index_t e = m_faceOffset[5]+INDEX2(k0,k1,m_NE0);
10075                            /* GENERATOR SNIP_PDEBC_SINGLE_5 TOP */
10076                            ///////////////
10077                            // process d //
10078                            ///////////////
10079                            if (add_EM_S) {
10080                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10081                                if (d.actsExpanded()) {
10082                                    const double d_0 = d_p[0];
10083                                    const double d_1 = d_p[1];
10084                                    const double d_2 = d_p[2];
10085                                    const double d_3 = d_p[3];
10086                                    const double tmp2_0 = d_1 + d_3;
10087                                    const double tmp0_0 = d_0 + d_1 + d_2 + d_3;
10088                                    const double tmp1_0 = d_0 + d_2;
10089                                    const double tmp3_0 = d_0 + d_1;
10090                                    const double tmp6_0 = d_0 + d_3;
10091                                    const double tmp4_0 = d_2 + d_3;
10092                                    const double tmp5_0 = d_1 + d_2;
10093                                    const double tmp1_1 = tmp1_0*w25;
10094                                    const double tmp11_1 = d_0*w27;
10095                                    const double tmp9_1 = tmp5_0*w26;
10096                                    const double tmp12_1 = tmp2_0*w25;
10097                                    const double tmp3_1 = tmp3_0*w25;
10098                                    const double tmp6_1 = tmp3_0*w24;
10099                                    const double tmp2_1 = tmp2_0*w24;
10100                                    const double tmp10_1 = d_3*w28;
10101                                    const double tmp16_1 = tmp6_0*w26;
10102                                    const double tmp17_1 = d_1*w28;
10103                                    const double tmp7_1 = d_0*w28;
10104                                    const double tmp14_1 = d_2*w28;
10105                                    const double tmp18_1 = d_2*w27;
10106                                    const double tmp15_1 = d_1*w27;
10107                                    const double tmp0_1 = tmp0_0*w26;
10108                                    const double tmp4_1 = tmp4_0*w24;
10109                                    const double tmp5_1 = tmp4_0*w25;
10110                                    const double tmp8_1 = d_3*w27;
10111                                    const double tmp13_1 = tmp1_0*w24;
10112                                    EM_S[INDEX2(4,7,8)]+=tmp0_1;
10113                                    EM_S[INDEX2(6,4,8)]+=tmp1_1 + tmp2_1;
10114                                    EM_S[INDEX2(5,4,8)]+=tmp3_1 + tmp4_1;
10115                                    EM_S[INDEX2(5,6,8)]+=tmp0_1;
10116                                    EM_S[INDEX2(6,7,8)]+=tmp5_1 + tmp6_1;
10117                                    EM_S[INDEX2(7,6,8)]+=tmp5_1 + tmp6_1;
10118                                    EM_S[INDEX2(4,4,8)]+=tmp7_1 + tmp8_1 + tmp9_1;
10119                                    EM_S[INDEX2(7,7,8)]+=tmp10_1 + tmp11_1 + tmp9_1;
10120                                    EM_S[INDEX2(5,7,8)]+=tmp12_1 + tmp13_1;
10121                                    EM_S[INDEX2(6,6,8)]+=tmp14_1 + tmp15_1 + tmp16_1;
10122                                    EM_S[INDEX2(4,5,8)]+=tmp3_1 + tmp4_1;
10123                                    EM_S[INDEX2(5,5,8)]+=tmp16_1 + tmp17_1 + tmp18_1;
10124                                    EM_S[INDEX2(7,5,8)]+=tmp12_1 + tmp13_1;
10125                                    EM_S[INDEX2(6,5,8)]+=tmp0_1;
10126                                    EM_S[INDEX2(4,6,8)]+=tmp1_1 + tmp2_1;
10127                                    EM_S[INDEX2(7,4,8)]+=tmp0_1;
10128                                } else { /* constant data */
10129                                    const double tmp2_1 = d_p[0]*w31;
10130                                    const double tmp0_1 = d_p[0]*w30;
10131                                    const double tmp1_1 = d_p[0]*w29;
10132                                    EM_S[INDEX2(4,7,8)]+=tmp0_1;
10133                                    EM_S[INDEX2(6,4,8)]+=tmp1_1;
10134                                    EM_S[INDEX2(5,4,8)]+=tmp1_1;
10135                                    EM_S[INDEX2(5,6,8)]+=tmp0_1;
10136                                    EM_S[INDEX2(6,7,8)]+=tmp1_1;
10137                                    EM_S[INDEX2(7,6,8)]+=tmp1_1;
10138                                    EM_S[INDEX2(4,4,8)]+=tmp2_1;
10139                                    EM_S[INDEX2(7,7,8)]+=tmp2_1;
10140                                    EM_S[INDEX2(5,7,8)]+=tmp1_1;
10141                                    EM_S[INDEX2(6,6,8)]+=tmp2_1;
10142                                    EM_S[INDEX2(4,5,8)]+=tmp1_1;
10143                                    EM_S[INDEX2(5,5,8)]+=tmp2_1;
10144                                    EM_S[INDEX2(7,5,8)]+=tmp1_1;
10145                                    EM_S[INDEX2(6,5,8)]+=tmp0_1;
10146                                    EM_S[INDEX2(4,6,8)]+=tmp1_1;
10147                                    EM_S[INDEX2(7,4,8)]+=tmp0_1;
10148                                }
10149                            }
10150                            ///////////////
10151                            // process y //
10152                            ///////////////
10153                            if (add_EM_F) {
10154                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10155                                if (y.actsExpanded()) {
10156                                    const double y_0 = y_p[0];
10157                                    const double y_1 = y_p[1];
10158                                    const double y_2 = y_p[2];
10159                                    const double y_3 = y_p[3];
10160                                    const double tmp0_0 = y_1 + y_2;
10161                                    const double tmp1_0 = y_0 + y_3;
10162                                    const double tmp7_1 = w34*y_1;
10163                                    const double tmp3_1 = w32*y_1;
10164                                    const double tmp8_1 = w32*y_3;
10165                                    const double tmp0_1 = w32*y_0;
10166                                    const double tmp2_1 = w34*y_3;
10167                                    const double tmp9_1 = w34*y_0;
10168                                    const double tmp6_1 = w32*y_2;
10169                                    const double tmp5_1 = w34*y_2;
10170                                    const double tmp1_1 = tmp0_0*w33;
10171                                    const double tmp4_1 = tmp1_0*w33;
10172                                    EM_F[4]+=tmp0_1 + tmp1_1 + tmp2_1;
10173                                    EM_F[5]+=tmp3_1 + tmp4_1 + tmp5_1;
10174                                    EM_F[6]+=tmp4_1 + tmp6_1 + tmp7_1;
10175                                    EM_F[7]+=tmp1_1 + tmp8_1 + tmp9_1;
10176                                } else { /* constant data */
10177                                    const double tmp0_1 = w35*y_p[0];
10178                                    EM_F[4]+=tmp0_1;
10179                                    EM_F[5]+=tmp0_1;
10180                                    EM_F[6]+=tmp0_1;
10181                                    EM_F[7]+=tmp0_1;
10182                                }
10183                            }
10184                            /* GENERATOR SNIP_PDEBC_SINGLE_5 BOTTOM */
10185                            const index_t firstNode=m_N0*m_N1*(m_N2-2)+m_N0*k1+k0;
10186                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10187                                    add_EM_F, firstNode);
10188                        } // k0 loop
10189                    } // k1 loop
10190                } // colouring
10191            } // face 5
10192        } // end of parallel region
10193  }  }
10194    
10195  //protected  //protected
10196  void Brick::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,  void Brick::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
10197        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
10198  {  {
10199        const double h0 = m_l0/m_gNE0;
10200        const double h1 = m_l1/m_gNE1;
10201        const double h2 = m_l2/m_gNE2;
10202        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE TOP */
10203        const double w0 = 0.0625*h1*h2;
10204        const double w1 = 0.25*h1*h2;
10205        const double w2 = 0.0625*h0*h2;
10206        const double w3 = 0.25*h0*h2;
10207        const double w4 = 0.0625*h0*h1;
10208        const double w5 = 0.25*h0*h1;
10209        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE BOTTOM */
10210        const bool add_EM_S=!d.isEmpty();
10211        const bool add_EM_F=!y.isEmpty();
10212        rhs.requireWrite();
10213    #pragma omp parallel
10214        {
10215            if (m_faceOffset[0] > -1) {
10216                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10217    #pragma omp for nowait
10218                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10219                        for (index_t k1=0; k1<m_NE1; ++k1) {
10220                            vector<double> EM_S(8*8, 0);
10221                            vector<double> EM_F(8, 0);
10222                            const index_t e = INDEX2(k1,k2,m_NE1);
10223                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 TOP */
10224                            ///////////////
10225                            // process d //
10226                            ///////////////
10227                            if (add_EM_S) {
10228                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10229                                const double d_0 = d_p[0];
10230                                const double tmp0_1 = d_0*w0;
10231                                EM_S[INDEX2(6,4,8)]+=tmp0_1;
10232                                EM_S[INDEX2(2,6,8)]+=tmp0_1;
10233                                EM_S[INDEX2(4,0,8)]+=tmp0_1;
10234                                EM_S[INDEX2(2,0,8)]+=tmp0_1;
10235                                EM_S[INDEX2(4,4,8)]+=tmp0_1;
10236                                EM_S[INDEX2(2,2,8)]+=tmp0_1;
10237                                EM_S[INDEX2(0,0,8)]+=tmp0_1;
10238                                EM_S[INDEX2(6,6,8)]+=tmp0_1;
10239                                EM_S[INDEX2(0,4,8)]+=tmp0_1;
10240                                EM_S[INDEX2(6,0,8)]+=tmp0_1;
10241                                EM_S[INDEX2(4,2,8)]+=tmp0_1;
10242                                EM_S[INDEX2(4,6,8)]+=tmp0_1;
10243                                EM_S[INDEX2(0,2,8)]+=tmp0_1;
10244                                EM_S[INDEX2(0,6,8)]+=tmp0_1;
10245                                EM_S[INDEX2(6,2,8)]+=tmp0_1;
10246                                EM_S[INDEX2(2,4,8)]+=tmp0_1;
10247                            }
10248                            ///////////////
10249                            // process y //
10250                            ///////////////
10251                            if (add_EM_F) {
10252                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10253                                const double y_0 = y_p[0];
10254                                const double tmp0_1 = w1*y_0;
10255                                EM_F[0]+=tmp0_1;
10256                                EM_F[2]+=tmp0_1;
10257                                EM_F[4]+=tmp0_1;
10258                                EM_F[6]+=tmp0_1;
10259                            }
10260                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 BOTTOM */
10261                            const index_t firstNode=m_N0*m_N1*k2+m_N0*k1;
10262                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10263                                    add_EM_F, firstNode);
10264                        } // k1 loop
10265                    } // k2 loop
10266                } // colouring
10267            } // face 0
10268    
10269            if (m_faceOffset[1] > -1) {
10270                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10271    #pragma omp for nowait
10272                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10273                        for (index_t k1=0; k1<m_NE1; ++k1) {
10274                            vector<double> EM_S(8*8, 0);
10275                            vector<double> EM_F(8, 0);
10276                            const index_t e = m_faceOffset[1]+INDEX2(k1,k2,m_NE1);
10277                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 TOP */
10278                            ///////////////
10279                            // process d //
10280                            ///////////////
10281                            if (add_EM_S) {
10282                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10283                                const double d_0 = d_p[0];
10284                                const double tmp0_1 = d_0*w0;
10285                                EM_S[INDEX2(7,3,8)]+=tmp0_1;
10286                                EM_S[INDEX2(1,3,8)]+=tmp0_1;
10287                                EM_S[INDEX2(5,1,8)]+=tmp0_1;
10288                                EM_S[INDEX2(3,7,8)]+=tmp0_1;
10289                                EM_S[INDEX2(3,3,8)]+=tmp0_1;
10290                                EM_S[INDEX2(1,5,8)]+=tmp0_1;
10291                                EM_S[INDEX2(7,7,8)]+=tmp0_1;
10292                                EM_S[INDEX2(5,7,8)]+=tmp0_1;
10293                                EM_S[INDEX2(5,3,8)]+=tmp0_1;
10294                                EM_S[INDEX2(1,1,8)]+=tmp0_1;
10295                                EM_S[INDEX2(7,1,8)]+=tmp0_1;
10296                                EM_S[INDEX2(5,5,8)]+=tmp0_1;
10297                                EM_S[INDEX2(7,5,8)]+=tmp0_1;
10298                                EM_S[INDEX2(3,5,8)]+=tmp0_1;
10299                                EM_S[INDEX2(3,1,8)]+=tmp0_1;
10300                                EM_S[INDEX2(1,7,8)]+=tmp0_1;
10301                            }
10302                            ///////////////
10303                            // process y //
10304                            ///////////////
10305                            if (add_EM_F) {
10306                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10307                                const double y_0 = y_p[0];
10308                                const double tmp0_1 = w1*y_0;
10309                                EM_F[1]+=tmp0_1;
10310                                EM_F[3]+=tmp0_1;
10311                                EM_F[5]+=tmp0_1;
10312                                EM_F[7]+=tmp0_1;
10313                            }
10314                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 BOTTOM */
10315                            const index_t firstNode=m_N0*m_N1*k2+m_N0*(k1+1)-2;
10316                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10317                                    add_EM_F, firstNode);
10318                        } // k1 loop
10319                    } // k2 loop
10320                } // colouring
10321            } // face 1
10322    
10323            if (m_faceOffset[2] > -1) {
10324                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10325    #pragma omp for nowait
10326                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10327                        for (index_t k0=0; k0<m_NE0; ++k0) {
10328                            vector<double> EM_S(8*8, 0);
10329                            vector<double> EM_F(8, 0);
10330                            const index_t e = m_faceOffset[2]+INDEX2(k0,k2,m_NE0);
10331                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 TOP */
10332                            ///////////////
10333                            // process d //
10334                            ///////////////
10335                            if (add_EM_S) {
10336                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10337                                const double d_0 = d_p[0];
10338                                const double tmp0_1 = d_0*w2;
10339                                EM_S[INDEX2(5,4,8)]+=tmp0_1;
10340                                EM_S[INDEX2(5,1,8)]+=tmp0_1;
10341                                EM_S[INDEX2(4,0,8)]+=tmp0_1;
10342                                EM_S[INDEX2(4,4,8)]+=tmp0_1;
10343                                EM_S[INDEX2(1,5,8)]+=tmp0_1;
10344                                EM_S[INDEX2(4,1,8)]+=tmp0_1;
10345                                EM_S[INDEX2(1,1,8)]+=tmp0_1;
10346                                EM_S[INDEX2(0,0,8)]+=tmp0_1;
10347                                EM_S[INDEX2(5,0,8)]+=tmp0_1;
10348                                EM_S[INDEX2(4,5,8)]+=tmp0_1;
10349                                EM_S[INDEX2(0,4,8)]+=tmp0_1;
10350                                EM_S[INDEX2(5,5,8)]+=tmp0_1;
10351                                EM_S[INDEX2(1,4,8)]+=tmp0_1;
10352                                EM_S[INDEX2(1,0,8)]+=tmp0_1;
10353                                EM_S[INDEX2(0,1,8)]+=tmp0_1;
10354                                EM_S[INDEX2(0,5,8)]+=tmp0_1;
10355                            }
10356                            ///////////////
10357                            // process y //
10358                            ///////////////
10359                            if (add_EM_F) {
10360                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10361                                const double y_0 = y_p[0];
10362                                const double tmp0_1 = w3*y_0;
10363                                EM_F[0]+=tmp0_1;
10364                                EM_F[1]+=tmp0_1;
10365                                EM_F[4]+=tmp0_1;
10366                                EM_F[5]+=tmp0_1;
10367                            }
10368                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 BOTTOM */
10369                            const index_t firstNode=m_N0*m_N1*k2+k0;
10370                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10371                                    add_EM_F, firstNode);
10372                        } // k0 loop
10373                    } // k2 loop
10374                } // colouring
10375            } // face 2
10376    
10377            if (m_faceOffset[3] > -1) {
10378                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10379    #pragma omp for nowait
10380                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10381                        for (index_t k0=0; k0<m_NE0; ++k0) {
10382                            vector<double> EM_S(8*8, 0);
10383                            vector<double> EM_F(8, 0);
10384                            const index_t e = m_faceOffset[3]+INDEX2(k0,k2,m_NE0);
10385                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 TOP */
10386                            ///////////////
10387                            // process d //
10388                            ///////////////
10389                            if (add_EM_S) {
10390                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10391                                const double d_0 = d_p[0];
10392                                const double tmp0_1 = d_0*w2;
10393                                EM_S[INDEX2(7,3,8)]+=tmp0_1;
10394                                EM_S[INDEX2(2,6,8)]+=tmp0_1;
10395                                EM_S[INDEX2(3,7,8)]+=tmp0_1;
10396                                EM_S[INDEX2(7,2,8)]+=tmp0_1;
10397                                EM_S[INDEX2(6,7,8)]+=tmp0_1;
10398                                EM_S[INDEX2(3,3,8)]+=tmp0_1;
10399                                EM_S[INDEX2(7,6,8)]+=tmp0_1;
10400                                EM_S[INDEX2(6,3,8)]+=tmp0_1;
10401                                EM_S[INDEX2(3,6,8)]+=tmp0_1;
10402                                EM_S[INDEX2(2,2,8)]+=tmp0_1;
10403                                EM_S[INDEX2(7,7,8)]+=tmp0_1;
10404                                EM_S[INDEX2(2,7,8)]+=tmp0_1;
10405                                EM_S[INDEX2(3,2,8)]+=tmp0_1;
10406                                EM_S[INDEX2(6,6,8)]+=tmp0_1;
10407                                EM_S[INDEX2(2,3,8)]+=tmp0_1;
10408                                EM_S[INDEX2(6,2,8)]+=tmp0_1;
10409                            }
10410                            ///////////////
10411                            // process y //
10412                            ///////////////
10413                            if (add_EM_F) {
10414                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10415                                const double y_0 = y_p[0];
10416                                const double tmp0_1 = w3*y_0;
10417                                EM_F[2]+=tmp0_1;
10418                                EM_F[3]+=tmp0_1;
10419                                EM_F[6]+=tmp0_1;
10420                                EM_F[7]+=tmp0_1;
10421                            }
10422                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 BOTTOM */
10423                            const index_t firstNode=m_N0*m_N1*k2+m_N0*(m_N1-2)+k0;
10424                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10425                                    add_EM_F, firstNode);
10426                        } // k0 loop
10427                    } // k2 loop
10428                } // colouring
10429            } // face 3
10430    
10431            if (m_faceOffset[4] > -1) {
10432                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
10433    #pragma omp for nowait
10434                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
10435                        for (index_t k0=0; k0<m_NE0; ++k0) {
10436                            vector<double> EM_S(8*8, 0);
10437                            vector<double> EM_F(8, 0);
10438                            const index_t e = m_faceOffset[4]+INDEX2(k0,k1,m_NE0);
10439                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_4 TOP */
10440                            ///////////////
10441                            // process d //
10442                            ///////////////
10443                            if (add_EM_S) {
10444                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10445                                const double d_0 = d_p[0];
10446                                const double tmp0_1 = d_0*w4;
10447                                EM_S[INDEX2(1,3,8)]+=tmp0_1;
10448                                EM_S[INDEX2(3,0,8)]+=tmp0_1;
10449                                EM_S[INDEX2(0,3,8)]+=tmp0_1;
10450                                EM_S[INDEX2(1,2,8)]+=tmp0_1;
10451                                EM_S[INDEX2(3,3,8)]+=tmp0_1;
10452                                EM_S[INDEX2(2,0,8)]+=tmp0_1;
10453                                EM_S[INDEX2(2,2,8)]+=tmp0_1;
10454                                EM_S[INDEX2(1,1,8)]+=tmp0_1;
10455                                EM_S[INDEX2(3,2,8)]+=tmp0_1;
10456                                EM_S[INDEX2(0,0,8)]+=tmp0_1;
10457                                EM_S[INDEX2(2,3,8)]+=tmp0_1;
10458                                EM_S[INDEX2(2,1,8)]+=tmp0_1;
10459                                EM_S[INDEX2(1,0,8)]+=tmp0_1;
10460                                EM_S[INDEX2(0,1,8)]+=tmp0_1;
10461                                EM_S[INDEX2(3,1,8)]+=tmp0_1;
10462                                EM_S[INDEX2(0,2,8)]+=tmp0_1;
10463                            }
10464                            ///////////////
10465                            // process y //
10466                            ///////////////
10467                            if (add_EM_F) {
10468                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10469                                const double y_0 = y_p[0];
10470                                const double tmp0_1 = w5*y_0;
10471                                EM_F[0]+=tmp0_1;
10472                                EM_F[1]+=tmp0_1;
10473                                EM_F[2]+=tmp0_1;
10474                                EM_F[3]+=tmp0_1;
10475                            }
10476                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_4 BOTTOM */
10477                            const index_t firstNode=m_N0*k1+k0;
10478                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10479                                    add_EM_F, firstNode);
10480                        } // k0 loop
10481                    } // k1 loop
10482                } // colouring
10483            } // face 4
10484    
10485            if (m_faceOffset[5] > -1) {
10486                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
10487    #pragma omp for nowait
10488                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
10489                        for (index_t k0=0; k0<m_NE0; ++k0) {
10490                            vector<double> EM_S(8*8, 0);
10491                            vector<double> EM_F(8, 0);
10492                            const index_t e = m_faceOffset[5]+INDEX2(k0,k1,m_NE0);
10493                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_5 TOP */
10494                            ///////////////
10495                            // process d //
10496                            ///////////////
10497                            if (add_EM_S) {
10498                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10499                                const double d_0 = d_p[0];
10500                                const double tmp0_1 = d_0*w4;
10501                                EM_S[INDEX2(4,7,8)]+=tmp0_1;
10502                                EM_S[INDEX2(6,4,8)]+=tmp0_1;
10503                                EM_S[INDEX2(5,4,8)]+=tmp0_1;
10504                                EM_S[INDEX2(5,6,8)]+=tmp0_1;
10505                                EM_S[INDEX2(6,7,8)]+=tmp0_1;
10506                                EM_S[INDEX2(7,6,8)]+=tmp0_1;
10507                                EM_S[INDEX2(4,4,8)]+=tmp0_1;
10508                                EM_S[INDEX2(7,7,8)]+=tmp0_1;
10509                                EM_S[INDEX2(5,7,8)]+=tmp0_1;
10510                                EM_S[INDEX2(6,6,8)]+=tmp0_1;
10511                                EM_S[INDEX2(4,5,8)]+=tmp0_1;
10512                                EM_S[INDEX2(5,5,8)]+=tmp0_1;
10513                                EM_S[INDEX2(7,5,8)]+=tmp0_1;
10514                                EM_S[INDEX2(6,5,8)]+=tmp0_1;
10515                                EM_S[INDEX2(4,6,8)]+=tmp0_1;
10516                                EM_S[INDEX2(7,4,8)]+=tmp0_1;
10517                            }
10518                            ///////////////
10519                            // process y //
10520                            ///////////////
10521                            if (add_EM_F) {
10522                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10523                                const double y_0 = y_p[0];
10524                                const double tmp0_1 = w5*y_0;
10525                                EM_F[4]+=tmp0_1;
10526                                EM_F[5]+=tmp0_1;
10527                                EM_F[6]+=tmp0_1;
10528                                EM_F[7]+=tmp0_1;
10529                            }
10530                            /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_5 BOTTOM */
10531                            const index_t firstNode=m_N0*m_N1*(m_N2-2)+m_N0*k1+k0;
10532                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10533                                    add_EM_F, firstNode);
10534                        } // k0 loop
10535                    } // k1 loop
10536                } // colouring
10537            } // face 5
10538        } // end of parallel region
10539  }  }
10540    
10541  //protected  //protected
10542  void Brick::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,  void Brick::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
10543        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
10544  {  {
10545        const double h0 = m_l0/m_gNE0;
10546        const double h1 = m_l1/m_gNE1;
10547        const double h2 = m_l2/m_gNE2;
10548        dim_t numEq, numComp;
10549        if (!mat)
10550            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
10551        else {
10552            numEq=mat->logical_row_block_size;
10553            numComp=mat->logical_col_block_size;
10554        }
10555        /* GENERATOR SNIP_PDEBC_SYSTEM_PRE TOP */
10556        const double w0 = 0.0018607582807716854616*h1*h2;
10557        const double w1 = 0.025917019497006092316*h1*h2;
10558        const double w10 = 0.01116454968463011277*h1*h2;
10559        const double w11 = 0.25*h1*h2;
10560        const double w12 = 0.0018607582807716854616*h0*h2;
10561        const double w13 = 0.025917019497006092316*h0*h2;
10562        const double w14 = 0.0069444444444444444444*h0*h2;
10563        const double w15 = 0.00049858867864229740201*h0*h2;
10564        const double w16 = 0.09672363354357992482*h0*h2;
10565        const double w17 = 0.055555555555555555556*h0*h2;
10566        const double w18 = 0.11111111111111111111*h0*h2;
10567        const double w19 = 0.027777777777777777778*h0*h2;
10568        const double w2 = 0.0069444444444444444444*h1*h2;
10569        const double w20 = 0.1555021169820365539*h0*h2;
10570        const double w21 = 0.041666666666666666667*h0*h2;
10571        const double w22 = 0.01116454968463011277*h0*h2;
10572        const double w23 = 0.25*h0*h2;
10573        const double w24 = 0.0018607582807716854616*h0*h1;
10574        const double w25 = 0.025917019497006092316*h0*h1;
10575        const double w26 = 0.0069444444444444444444*h0*h1;
10576        const double w27 = 0.00049858867864229740201*h0*h1;
10577        const double w28 = 0.09672363354357992482*h0*h1;
10578        const double w29 = 0.055555555555555555556*h0*h1;
10579        const double w3 = 0.00049858867864229740201*h1*h2;
10580        const double w30 = 0.027777777777777777778*h0*h1;
10581        const double w31 = 0.11111111111111111111*h0*h1;
10582        const double w32 = 0.1555021169820365539*h0*h1;
10583        const double w33 = 0.041666666666666666667*h0*h1;
10584        const double w34 = 0.01116454968463011277*h0*h1;
10585        const double w35 = 0.25*h0*h1;
10586        const double w4 = 0.09672363354357992482*h1*h2;
10587        const double w5 = 0.055555555555555555556*h1*h2;
10588        const double w6 = 0.11111111111111111111*h1*h2;
10589        const double w7 = 0.027777777777777777778*h1*h2;
10590        const double w8 = 0.1555021169820365539*h1*h2;
10591        const double w9 = 0.041666666666666666667*h1*h2;
10592        /* GENERATOR SNIP_PDEBC_SYSTEM_PRE BOTTOM */
10593        const bool add_EM_S=!d.isEmpty();
10594        const bool add_EM_F=!y.isEmpty();
10595        rhs.requireWrite();
10596    #pragma omp parallel
10597        {
10598            if (m_faceOffset[0] > -1) {
10599                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10600    #pragma omp for nowait
10601                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10602                        for (index_t k1=0; k1<m_NE1; ++k1) {
10603                            vector<double> EM_S(8*8*numEq*numComp, 0);
10604                            vector<double> EM_F(8*numEq, 0);
10605                            const index_t e = INDEX2(k1,k2,m_NE1);
10606                            /* GENERATOR SNIP_PDEBC_SYSTEM_0 TOP */
10607                            ///////////////
10608                            // process d //
10609                            ///////////////
10610                            if (add_EM_S) {
10611                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10612                                if (d.actsExpanded()) {
10613                                    for (index_t k=0; k<numEq; k++) {
10614                                        for (index_t m=0; m<numComp; m++) {
10615                                            const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
10616                                            const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
10617                                            const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];
10618                                            const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];
10619                                            const double tmp3_0 = d_1 + d_3;
10620                                            const double tmp6_0 = d_0 + d_1 + d_2 + d_3;
10621                                            const double tmp2_0 = d_0 + d_2;
10622                                            const double tmp0_0 = d_0 + d_1;
10623                                            const double tmp4_0 = d_0 + d_3;
10624                                            const double tmp1_0 = d_2 + d_3;
10625                                            const double tmp5_0 = d_1 + d_2;
10626                                            const double tmp14_1 = d_0*w4;
10627                                            const double tmp17_1 = d_3*w4;
10628                                            const double tmp12_1 = d_1*w4;
10629                                            const double tmp9_1 = d_2*w4;
10630                                            const double tmp4_1 = tmp3_0*w0;
10631                                            const double tmp3_1 = tmp3_0*w1;
10632                                            const double tmp7_1 = tmp0_0*w1;
10633                                            const double tmp8_1 = d_1*w3;
10634                                            const double tmp16_1 = d_0*w3;
10635                                            const double tmp18_1 = tmp6_0*w2;
10636                                            const double tmp1_1 = tmp1_0*w1;
10637                                            const double tmp15_1 = tmp5_0*w2;
10638                                            const double tmp6_1 = tmp1_0*w0;
10639                                            const double tmp2_1 = tmp2_0*w0;
10640                                            const double tmp13_1 = d_3*w3;
10641                                            const double tmp5_1 = tmp2_0*w1;
10642                                            const double tmp10_1 = tmp4_0*w2;
10643                                            const double tmp11_1 = d_2*w3;
10644                                            const double tmp0_1 = tmp0_0*w0;
10645                                            EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
10646                                            EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp2_1 + tmp3_1;
10647                                            EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp4_1 + tmp5_1;
10648                                            EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp6_1 + tmp7_1;
10649                                            EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp10_1 + tmp8_1 + tmp9_1;
10650                                            EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp12_1;
10651                                            EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp13_1 + tmp14_1 + tmp15_1;
10652                                            EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp15_1 + tmp16_1 + tmp17_1;
10653                                            EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp4_1 + tmp5_1;
10654                                            EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp18_1;
10655                                            EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp18_1;
10656                                            EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
10657                                            EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp6_1 + tmp7_1;
10658                                            EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp18_1;
10659                                            EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp2_1 + tmp3_1;
10660                                            EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp18_1;
10661                                        }
10662                                     }
10663                                } else { /* constant data */
10664                                    for (index_t k=0; k<numEq; k++) {
10665                                        for (index_t m=0; m<numComp; m++) {
10666                                            const double d_0 = d_p[INDEX2(k, m, numEq)];
10667                                            const double tmp0_1 = d_0*w5;
10668                                            const double tmp2_1 = d_0*w7;
10669                                            const double tmp1_1 = d_0*w6;
10670                                            EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0_1;
10671                                            EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1;
10672                                            EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0_1;
10673                                            EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0_1;
10674                                            EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp1_1;
10675                                            EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp1_1;
10676                                            EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp1_1;
10677                                            EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp1_1;
10678                                            EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0_1;
10679                                            EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp2_1;
10680                                            EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp2_1;
10681                                            EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1;
10682                                            EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1;
10683                                            EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp2_1;
10684                                            EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1;
10685                                            EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp2_1;
10686                                        }
10687                                    }
10688                                }
10689                            }
10690                            ///////////////
10691                            // process y //
10692                            ///////////////
10693                            if (add_EM_F) {
10694                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10695                                if (y.actsExpanded()) {
10696                                    for (index_t k=0; k<numEq; k++) {
10697                                        const double y_0 = y_p[INDEX2(k, 0, numEq)];
10698                                        const double y_1 = y_p[INDEX2(k, 1, numEq)];
10699                                        const double y_2 = y_p[INDEX2(k, 2, numEq)];
10700                                        const double y_3 = y_p[INDEX2(k, 3, numEq)];
10701                                        const double tmp0_0 = y_1 + y_2;
10702                                        const double tmp1_0 = y_0 + y_3;
10703                                        const double tmp2_1 = w10*y_3;
10704                                        const double tmp8_1 = w10*y_0;
10705                                        const double tmp5_1 = w10*y_2;
10706                                        const double tmp3_1 = w8*y_1;
10707                                        const double tmp9_1 = w8*y_3;
10708                                        const double tmp0_1 = w8*y_0;
10709                                        const double tmp1_1 = tmp0_0*w9;
10710                                        const double tmp7_1 = w8*y_2;
10711                                        const double tmp4_1 = tmp1_0*w9;
10712                                        const double tmp6_1 = w10*y_1;
10713                                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
10714                                        EM_F[INDEX2(k,2,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
10715                                        EM_F[INDEX2(k,4,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
10716                                        EM_F[INDEX2(k,6,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;
10717                                    }
10718                                } else { /* constant data */
10719                                    for (index_t k=0; k<numEq; k++) {
10720                                        const double y_0 = y_p[k];
10721                                        const double tmp0_1 = w11*y_0;
10722                                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
10723                                        EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
10724                                        EM_F[INDEX2(k,4,numEq)]+=tmp0_1;
10725                                        EM_F[INDEX2(k,6,numEq)]+=tmp0_1;
10726                                    }
10727                                }
10728                            }
10729                            /* GENERATOR SNIP_PDEBC_SYSTEM_0 BOTTOM */
10730                            const index_t firstNode=m_N0*m_N1*k2+m_N0*k1;
10731                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10732                                    add_EM_F, firstNode, numEq, numComp);
10733                        } // k1 loop
10734                    } // k2 loop
10735                } // colouring
10736            } // face 0
10737    
10738            if (m_faceOffset[1] > -1) {
10739                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10740    #pragma omp for nowait
10741                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10742                        for (index_t k1=0; k1<m_NE1; ++k1) {
10743                            vector<double> EM_S(8*8*numEq*numComp, 0);
10744                            vector<double> EM_F(8*numEq, 0);
10745                            const index_t e = m_faceOffset[1]+INDEX2(k1,k2,m_NE1);
10746                            /* GENERATOR SNIP_PDEBC_SYSTEM_1 TOP */
10747                            ///////////////
10748                            // process d //
10749                            ///////////////
10750                            if (add_EM_S) {
10751                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10752                                if (d.actsExpanded()) {
10753                                    for (index_t k=0; k<numEq; k++) {
10754                                        for (index_t m=0; m<numComp; m++) {
10755                                            const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
10756                                            const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
10757                                            const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];
10758                                            const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];
10759                                            const double tmp1_0 = d_1 + d_3;
10760                                            const double tmp6_0 = d_0 + d_1 + d_2 + d_3;
10761                                            const double tmp0_0 = d_0 + d_2;
10762                                            const double tmp3_0 = d_0 + d_1;
10763                                            const double tmp4_0 = d_0 + d_3;
10764                                            const double tmp2_0 = d_2 + d_3;
10765                                            const double tmp5_0 = d_1 + d_2;
10766                                            const double tmp10_1 = d_3*w4;
10767                                            const double tmp13_1 = tmp2_0*w1;
10768                                            const double tmp16_1 = d_0*w4;
10769                                            const double tmp18_1 = d_2*w4;
10770                                            const double tmp12_1 = tmp3_0*w0;
10771                                            const double tmp3_1 = tmp3_0*w1;
10772                                            const double tmp5_1 = tmp0_0*w1;
10773                                            const double tmp17_1 = d_1*w3;
10774                                            const double tmp9_1 = d_0*w3;
10775                                            const double tmp14_1 = tmp6_0*w2;
10776                                            const double tmp1_1 = tmp1_0*w1;
10777                                            const double tmp11_1 = tmp5_0*w2;
10778                                            const double tmp4_1 = tmp1_0*w0;
10779                                            const double tmp2_1 = tmp2_0*w0;
10780                                            const double tmp15_1 = d_3*w3;
10781                                            const double tmp7_1 = d_1*w4;
10782                                            const double tmp8_1 = tmp4_0*w2;
10783                                            const double tmp6_1 = d_2*w3;
10784                                            const double tmp0_1 = tmp0_0*w0;
10785                                            EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
10786                                            EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp2_1 + tmp3_1;
10787                                            EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp4_1 + tmp5_1;
10788                                            EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
10789                                            EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp6_1 + tmp7_1 + tmp8_1;
10790                                            EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp4_1 + tmp5_1;
10791                                            EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp9_1;
10792                                            EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp12_1 + tmp13_1;
10793                                            EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp14_1;
10794                                            EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp11_1 + tmp15_1 + tmp16_1;
10795                                            EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp14_1;
10796                                            EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp17_1 + tmp18_1 + tmp8_1;
10797                                            EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp12_1 + tmp13_1;
10798                                            EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp14_1;
10799                                            EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp2_1 + tmp3_1;
10800                                            EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp14_1;
10801                                        }
10802                                     }
10803                                } else { /* constant data */
10804                                    for (index_t k=0; k<numEq; k++) {
10805                                        for (index_t m=0; m<numComp; m++) {
10806                                            const double d_0 = d_p[INDEX2(k, m, numEq)];
10807                                            const double tmp0_1 = d_0*w5;
10808                                            const double tmp2_1 = d_0*w7;
10809                                            const double tmp1_1 = d_0*w6;
10810                                            EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1;
10811                                            EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1;
10812                                            EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0_1;
10813                                            EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1;
10814                                            EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp1_1;
10815                                            EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0_1;
10816                                            EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp1_1;
10817                                            EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0_1;
10818                                            EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp2_1;
10819                                            EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp1_1;
10820                                            EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp2_1;
10821                                            EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp1_1;
10822                                            EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0_1;
10823                                            EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp2_1;
10824                                            EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1;
10825                                            EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp2_1;
10826                                        }
10827                                    }
10828                                }
10829                            }
10830                            ///////////////
10831                            // process y //
10832                            ///////////////
10833                            if (add_EM_F) {
10834                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10835                                if (y.actsExpanded()) {
10836                                    for (index_t k=0; k<numEq; k++) {
10837                                        const double y_0 = y_p[INDEX2(k, 0, numEq)];
10838                                        const double y_1 = y_p[INDEX2(k, 1, numEq)];
10839                                        const double y_2 = y_p[INDEX2(k, 2, numEq)];
10840                                        const double y_3 = y_p[INDEX2(k, 3, numEq)];
10841                                        const double tmp0_0 = y_1 + y_2;
10842                                        const double tmp1_0 = y_0 + y_3;
10843                                        const double tmp2_1 = w10*y_3;
10844                                        const double tmp8_1 = w10*y_0;
10845                                        const double tmp5_1 = w10*y_2;
10846                                        const double tmp3_1 = w8*y_1;
10847                                        const double tmp9_1 = w8*y_3;
10848                                        const double tmp0_1 = w8*y_0;
10849                                        const double tmp1_1 = tmp0_0*w9;
10850                                        const double tmp7_1 = w8*y_2;
10851                                        const double tmp4_1 = tmp1_0*w9;
10852                                        const double tmp6_1 = w10*y_1;
10853                                        EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
10854                                        EM_F[INDEX2(k,3,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
10855                                        EM_F[INDEX2(k,5,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
10856                                        EM_F[INDEX2(k,7,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;
10857                                    }
10858                                } else { /* constant data */
10859                                    for (index_t k=0; k<numEq; k++) {
10860                                        const double y_0 = y_p[k];
10861                                        const double tmp0_1 = w11*y_0;
10862                                        EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
10863                                        EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
10864                                        EM_F[INDEX2(k,5,numEq)]+=tmp0_1;
10865                                        EM_F[INDEX2(k,7,numEq)]+=tmp0_1;
10866                                    }
10867                                }
10868                            }
10869                            /* GENERATOR SNIP_PDEBC_SYSTEM_1 BOTTOM */
10870                            const index_t firstNode=m_N0*m_N1*k2+m_N0*(k1+1)-2;
10871                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
10872                                    add_EM_F, firstNode, numEq, numComp);
10873                        } // k1 loop
10874                    } // k2 loop
10875                } // colouring
10876            } // face 1
10877    
10878            if (m_faceOffset[2] > -1) {
10879                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10880    #pragma omp for nowait
10881                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10882                        for (index_t k0=0; k0<m_NE0; ++k0) {
10883                            vector<double> EM_S(8*8*numEq*numComp, 0);
10884                            vector<double> EM_F(8*numEq, 0);
10885                            const index_t e = m_faceOffset[2]+INDEX2(k0,k2,m_NE0);
10886                            /* GENERATOR SNIP_PDEBC_SYSTEM_2 TOP */
10887                            ///////////////
10888                            // process d //
10889                            ///////////////
10890                            if (add_EM_S) {
10891                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
10892                                if (d.actsExpanded()) {
10893                                    for (index_t k=0; k<numEq; k++) {
10894                                        for (index_t m=0; m<numComp; m++) {
10895                                            const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
10896                                            const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
10897                                            const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];
10898                                            const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];
10899                                            const double tmp2_0 = d_1 + d_3;
10900                                            const double tmp5_0 = d_0 + d_1 + d_2 + d_3;
10901                                            const double tmp3_0 = d_0 + d_2;
10902                                            const double tmp1_0 = d_0 + d_1;
10903                                            const double tmp4_0 = d_0 + d_3;
10904                                            const double tmp0_0 = d_2 + d_3;
10905                                            const double tmp6_0 = d_1 + d_2;
10906                                            const double tmp2_1 = tmp2_0*w13;
10907                                            const double tmp14_1 = d_3*w15;
10908                                            const double tmp0_1 = tmp0_0*w13;
10909                                            const double tmp3_1 = tmp3_0*w12;
10910                                            const double tmp17_1 = tmp1_0*w13;
10911                                            const double tmp18_1 = tmp0_0*w12;
10912                                            const double tmp8_1 = d_1*w15;
10913                                            const double tmp16_1 = d_0*w15;
10914                                            const double tmp11_1 = d_2*w15;
10915                                            const double tmp5_1 = tmp2_0*w12;
10916                                            const double tmp15_1 = d_3*w16;
10917                                            const double tmp13_1 = tmp6_0*w14;
10918                                            const double tmp1_1 = tmp1_0*w12;
10919                                            const double tmp7_1 = tmp4_0*w14;
10920                                            const double tmp12_1 = d_0*w16;
10921                                            const double tmp10_1 = d_1*w16;
10922                                            const double tmp6_1 = d_2*w16;
10923                                            const double tmp9_1 = tmp5_0*w14;
10924                                            const double tmp4_1 = tmp3_0*w13;
10925                                            EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
10926                                            EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp2_1 + tmp3_1;
10927                                            EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp4_1 + tmp5_1;
10928                                            EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp6_1 + tmp7_1 + tmp8_1;
10929                                            EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1;
10930                                            EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp9_1;
10931                                            EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp7_1;
10932                                            EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp12_1 + tmp13_1 + tmp14_1;
10933                                            EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp9_1;
10934                                            EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
10935                                            EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp4_1 + tmp5_1;
10936                                            EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp13_1 + tmp15_1 + tmp16_1;
10937                                            EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp9_1;
10938                                            EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp17_1 + tmp18_1;
10939                                            EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp17_1 + tmp18_1;
10940                                            EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp9_1;
10941                                        }
10942                                     }
10943                                } else { /* constant data */
10944                                    for (index_t k=0; k<numEq; k++) {
10945                                        for (index_t m=0; m<numComp; m++) {
10946                                            const double d_0 = d_p[INDEX2(k, m, numEq)];
10947                                            const double tmp0_1 = d_0*w17;
10948                                            const double tmp2_1 = d_0*w19;
10949                                            const double tmp1_1 = d_0*w18;
10950                                            EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0_1;
10951                                            EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0_1;
10952                                            EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0_1;
10953                                            EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp1_1;
10954                                            EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0_1;
10955                                            EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp2_1;
10956                                            EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp1_1;
10957                                            EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp1_1;
10958                                            EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp2_1;
10959                                            EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0_1;
10960                                            EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0_1;
10961                                            EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp1_1;
10962                                            EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp2_1;
10963                                            EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0_1;
10964                                            EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0_1;
10965                                            EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp2_1;
10966                                        }
10967                                    }
10968                                }
10969                            }
10970                            ///////////////
10971                            // process y //
10972                            ///////////////
10973                            if (add_EM_F) {
10974                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
10975                                if (y.actsExpanded()) {
10976                                    for (index_t k=0; k<numEq; k++) {
10977                                        const double y_0 = y_p[INDEX2(k, 0, numEq)];
10978                                        const double y_1 = y_p[INDEX2(k, 1, numEq)];
10979                                        const double y_2 = y_p[INDEX2(k, 2, numEq)];
10980                                        const double y_3 = y_p[INDEX2(k, 3, numEq)];
10981                                        const double tmp0_0 = y_1 + y_2;
10982                                        const double tmp1_0 = y_0 + y_3;
10983                                        const double tmp0_1 = w22*y_3;
10984                                        const double tmp6_1 = w22*y_1;
10985                                        const double tmp3_1 = w22*y_2;
10986                                        const double tmp5_1 = w20*y_1;
10987                                        const double tmp9_1 = w20*y_3;
10988                                        const double tmp4_1 = tmp1_0*w21;
10989                                        const double tmp8_1 = w22*y_0;
10990                                        const double tmp2_1 = w20*y_0;
10991                                        const double tmp7_1 = w20*y_2;
10992                                        const double tmp1_1 = tmp0_0*w21;
10993                                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
10994                                        EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
10995                                        EM_F[INDEX2(k,4,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
10996                                        EM_F[INDEX2(k,5,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;
10997                                    }
10998                                } else { /* constant data */
10999                                    for (index_t k=0; k<numEq; k++) {
11000                                        const double y_0 = y_p[k];
11001                                        const double tmp0_1 = w23*y_0;
11002                                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
11003                                        EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
11004                                        EM_F[INDEX2(k,4,numEq)]+=tmp0_1;
11005                                        EM_F[INDEX2(k,5,numEq)]+=tmp0_1;
11006                                    }
11007                                }
11008                            }
11009                            /* GENERATOR SNIP_PDEBC_SYSTEM_2 BOTTOM */
11010                            const index_t firstNode=m_N0*m_N1*k2+k0;
11011                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11012                                    add_EM_F, firstNode, numEq, numComp);
11013                        } // k0 loop
11014                    } // k2 loop
11015                } // colouring
11016            } // face 2
11017    
11018            if (m_faceOffset[3] > -1) {
11019                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11020    #pragma omp for nowait
11021                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11022                        for (index_t k0=0; k0<m_NE0; ++k0) {
11023                            vector<double> EM_S(8*8*numEq*numComp, 0);
11024                            vector<double> EM_F(8*numEq, 0);
11025                            const index_t e = m_faceOffset[3]+INDEX2(k0,k2,m_NE0);
11026                            /* GENERATOR SNIP_PDEBC_SYSTEM_3 TOP */
11027                            ///////////////
11028                            // process d //
11029                            ///////////////
11030                            if (add_EM_S) {
11031                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
11032                                if (d.actsExpanded()) {
11033                                    for (index_t k=0; k<numEq; k++) {
11034                                        for (index_t m=0; m<numComp; m++) {
11035                                            const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
11036                                            const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
11037                                            const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];
11038                                            const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];
11039                                            const double tmp0_0 = d_1 + d_3;
11040                                            const double tmp2_0 = d_0 + d_1 + d_2 + d_3;
11041                                            const double tmp1_0 = d_0 + d_2;
11042                                            const double tmp4_0 = d_0 + d_1;
11043                                            const double tmp5_0 = d_0 + d_3;
11044                                            const double tmp3_0 = d_2 + d_3;
11045                                            const double tmp6_0 = d_1 + d_2;
11046                                            const double tmp15_1 = tmp4_0*w13;
11047                                            const double tmp10_1 = d_0*w16;
11048                                            const double tmp6_1 = tmp4_0*w12;
11049                                            const double tmp16_1 = tmp3_0*w12;
11050                                            const double tmp0_1 = tmp0_0*w13;
11051                                            const double tmp2_1 = tmp1_0*w13;
11052                                            const double tmp18_1 = d_1*w15;
11053                                            const double tmp14_1 = d_0*w15;
11054                                            const double tmp9_1 = d_2*w15;
11055                                            const double tmp4_1 = tmp2_0*w14;
11056                                            const double tmp13_1 = d_3*w16;
11057                                            const double tmp11_1 = tmp6_0*w14;
11058                                            const double tmp1_1 = tmp1_0*w12;
11059                                            const double tmp12_1 = d_3*w15;
11060                                            const double tmp3_1 = tmp0_0*w12;
11061                                            const double tmp7_1 = d_1*w16;
11062                                            const double tmp17_1 = d_2*w16;
11063                                            const double tmp8_1 = tmp5_0*w14;
11064                                            const double tmp5_1 = tmp3_0*w13;
11065                                            EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
11066                                            EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp2_1 + tmp3_1;
11067                                            EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
11068                                            EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp4_1;
11069                                            EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp5_1 + tmp6_1;
11070                                            EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp7_1 + tmp8_1 + tmp9_1;
11071                                            EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp5_1 + tmp6_1;
11072                                            EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp4_1;
11073                                            EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp4_1;
11074                                            EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp12_1;
11075                                            EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp11_1 + tmp13_1 + tmp14_1;
11076                                            EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp4_1;
11077                                            EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp15_1 + tmp16_1;
11078                                            EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp17_1 + tmp18_1 + tmp8_1;
11079                                            EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp15_1 + tmp16_1;
11080                                            EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp2_1 + tmp3_1;
11081                                        }
11082                                     }
11083                                } else { /* constant data */
11084                                    for (index_t k=0; k<numEq; k++) {
11085                                        for (index_t m=0; m<numComp; m++) {
11086                                            const double d_0 = d_p[INDEX2(k, m, numEq)];
11087                                            const double tmp0_1 = d_0*w17;
11088                                            const double tmp1_1 = d_0*w19;
11089                                            const double tmp2_1 = d_0*w18;
11090                                            EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1;
11091                                            EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1;
11092                                            EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1;
11093                                            EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp1_1;
11094                                            EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0_1;
11095                                            EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp2_1;
11096                                            EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0_1;
11097                                            EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp1_1;
11098                                            EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp1_1;
11099                                            EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp2_1;
11100                                            EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp2_1;
11101                                            EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp1_1;
11102                                            EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1;
11103                                            EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp2_1;
11104                                            EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1;
11105                                            EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1;
11106                                        }
11107                                    }
11108                                }
11109                            }
11110                            ///////////////
11111                            // process y //
11112                            ///////////////
11113                            if (add_EM_F) {
11114                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
11115                                if (y.actsExpanded()) {
11116                                    for (index_t k=0; k<numEq; k++) {
11117                                        const double y_0 = y_p[INDEX2(k, 0, numEq)];
11118                                        const double y_1 = y_p[INDEX2(k, 1, numEq)];
11119                                        const double y_2 = y_p[INDEX2(k, 2, numEq)];
11120                                        const double y_3 = y_p[INDEX2(k, 3, numEq)];
11121                                        const double tmp0_0 = y_1 + y_2;
11122                                        const double tmp1_0 = y_0 + y_3;
11123                                        const double tmp0_1 = w22*y_3;
11124                                        const double tmp6_1 = w22*y_1;
11125                                        const double tmp3_1 = w22*y_2;
11126                                        const double tmp5_1 = w20*y_1;
11127                                        const double tmp9_1 = w20*y_3;
11128                                        const double tmp4_1 = tmp1_0*w21;
11129                                        const double tmp8_1 = w22*y_0;
11130                                        const double tmp2_1 = w20*y_0;
11131                                        const double tmp7_1 = w20*y_2;
11132                                        const double tmp1_1 = tmp0_0*w21;
11133                                        EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
11134                                        EM_F[INDEX2(k,3,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
11135                                        EM_F[INDEX2(k,6,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
11136                                        EM_F[INDEX2(k,7,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;
11137                                    }
11138                                } else { /* constant data */
11139                                    for (index_t k=0; k<numEq; k++) {
11140                                        const double y_0 = y_p[k];
11141                                        const double tmp0_1 = w23*y_0;
11142                                        EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
11143                                        EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
11144                                        EM_F[INDEX2(k,6,numEq)]+=tmp0_1;
11145                                        EM_F[INDEX2(k,7,numEq)]+=tmp0_1;
11146                                    }
11147                                }
11148                            }
11149                            /* GENERATOR SNIP_PDEBC_SYSTEM_3 BOTTOM */
11150                            const index_t firstNode=m_N0*m_N1*k2+m_N0*(m_N1-2)+k0;
11151                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11152                                    add_EM_F, firstNode, numEq, numComp);
11153                        } // k0 loop
11154                    } // k2 loop
11155                } // colouring
11156            } // face 3
11157    
11158            if (m_faceOffset[4] > -1) {
11159                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
11160    #pragma omp for nowait
11161                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
11162                        for (index_t k0=0; k0<m_NE0; ++k0) {
11163                            vector<double> EM_S(8*8*numEq*numComp, 0);
11164                            vector<double> EM_F(8*numEq, 0);
11165                            const index_t e = m_faceOffset[4]+INDEX2(k0,k1,m_NE0);
11166                            /* GENERATOR SNIP_PDEBC_SYSTEM_4 TOP */
11167                            ///////////////
11168                            // process d //
11169                            ///////////////
11170                            if (add_EM_S) {
11171                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
11172                                if (d.actsExpanded()) {
11173                                    for (index_t k=0; k<numEq; k++) {
11174                                        for (index_t m=0; m<numComp; m++) {
11175                                            const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
11176                                            const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
11177                                            const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];
11178                                            const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];
11179                                            const double tmp0_0 = d_1 + d_3;
11180                                            const double tmp2_0 = d_0 + d_1 + d_2 + d_3;
11181                                            const double tmp1_0 = d_0 + d_2;
11182                                            const double tmp6_0 = d_0 + d_1;
11183                                            const double tmp4_0 = d_0 + d_3;
11184                                            const double tmp5_0 = d_2 + d_3;
11185                                            const double tmp3_0 = d_1 + d_2;
11186                                            const double tmp18_1 = tmp5_0*w24;
11187                                            const double tmp6_1 = tmp1_0*w25;
11188                                            const double tmp4_1 = d_0*w27;
11189                                            const double tmp12_1 = d_2*w27;
11190                                            const double tmp0_1 = tmp0_0*w25;
11191                                            const double tmp5_1 = tmp3_0*w26;
11192                                            const double tmp2_1 = tmp2_0*w26;
11193                                            const double tmp17_1 = tmp6_0*w25;
11194                                            const double tmp14_1 = tmp6_0*w24;
11195                                            const double tmp11_1 = d_1*w28;
11196                                            const double tmp9_1 = d_1*w27;
11197                                            const double tmp16_1 = d_3*w27;
11198                                            const double tmp8_1 = d_2*w28;
11199                                            const double tmp7_1 = tmp0_0*w24;
11200                                            const double tmp15_1 = d_0*w28;
11201                                            const double tmp13_1 = tmp5_0*w25;
11202                                            const double tmp3_1 = d_3*w28;
11203                                            const double tmp10_1 = tmp4_0*w26;
11204                                            const double tmp1_1 = tmp1_0*w24;
11205                                            EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
11206                                            EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp2_1;
11207                                            EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp2_1;
11208                                            EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp2_1;
11209                                            EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
11210                                            EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp6_1 + tmp7_1;
11211                                            EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp10_1 + tmp8_1 + tmp9_1;
11212                                            EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp12_1;
11213                                            EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp13_1 + tmp14_1;
11214                                            EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp15_1 + tmp16_1 + tmp5_1;
11215                                            EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp13_1 + tmp14_1;
11216                                            EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp2_1;
11217                                            EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp17_1 + tmp18_1;
11218                                            EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp17_1 + tmp18_1;
11219                                            EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1 + tmp1_1;
11220                                            EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp6_1 + tmp7_1;
11221                                        }
11222                                     }
11223                                } else { /* constant data */
11224                                    for (index_t k=0; k<numEq; k++) {
11225                                        for (index_t m=0; m<numComp; m++) {
11226                                            const double d_0 = d_p[INDEX2(k, m, numEq)];
11227                                            const double tmp2_1 = d_0*w31;
11228                                            const double tmp1_1 = d_0*w30;
11229                                            const double tmp0_1 = d_0*w29;
11230                                            EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1;
11231                                            EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp1_1;
11232                                            EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp1_1;
11233                                            EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp1_1;
11234                                            EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp2_1;
11235                                            EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0_1;
11236                                            EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp2_1;
11237                                            EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp2_1;
11238                                            EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1;
11239                                            EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp2_1;
11240                                            EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1;
11241                                            EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp1_1;
11242                                            EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0_1;
11243                                            EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0_1;
11244                                            EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1;
11245                                            EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1;
11246                                        }
11247                                    }
11248                                }
11249                            }
11250                            ///////////////
11251                            // process y //
11252                            ///////////////
11253                            if (add_EM_F) {
11254                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
11255                                if (y.actsExpanded()) {
11256                                    for (index_t k=0; k<numEq; k++) {
11257                                        const double y_0 = y_p[INDEX2(k, 0, numEq)];
11258                                        const double y_1 = y_p[INDEX2(k, 1, numEq)];
11259                                        const double y_2 = y_p[INDEX2(k, 2, numEq)];
11260                                        const double y_3 = y_p[INDEX2(k, 3, numEq)];
11261                                        const double tmp0_0 = y_1 + y_2;
11262                                        const double tmp1_0 = y_0 + y_3;
11263                                        const double tmp7_1 = w34*y_1;
11264                                        const double tmp3_1 = w32*y_1;
11265                                        const double tmp8_1 = w32*y_3;
11266                                        const double tmp0_1 = w32*y_0;
11267                                        const double tmp2_1 = w34*y_3;
11268                                        const double tmp9_1 = w34*y_0;
11269                                        const double tmp6_1 = w32*y_2;
11270                                        const double tmp5_1 = w34*y_2;
11271                                        const double tmp1_1 = tmp0_0*w33;
11272                                        const double tmp4_1 = tmp1_0*w33;
11273                                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
11274                                        EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
11275                                        EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
11276                                        EM_F[INDEX2(k,3,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;
11277                                    }
11278                                } else { /* constant data */
11279                                    for (index_t k=0; k<numEq; k++) {
11280                                        const double y_0 = y_p[k];
11281                                        const double tmp0_1 = w35*y_0;
11282                                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
11283                                        EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
11284                                        EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
11285                                        EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
11286                                    }
11287                                }
11288                            }
11289                            /* GENERATOR SNIP_PDEBC_SYSTEM_4 BOTTOM */
11290                            const index_t firstNode=m_N0*k1+k0;
11291                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11292                                    add_EM_F, firstNode, numEq, numComp);
11293                        } // k0 loop
11294                    } // k1 loop
11295                } // colouring
11296            } // face 4
11297    
11298            if (m_faceOffset[5] > -1) {
11299                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
11300    #pragma omp for nowait
11301                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
11302                        for (index_t k0=0; k0<m_NE0; ++k0) {
11303                            vector<double> EM_S(8*8*numEq*numComp, 0);
11304                            vector<double> EM_F(8*numEq, 0);
11305                            const index_t e = m_faceOffset[5]+INDEX2(k0,k1,m_NE0);
11306                            /* GENERATOR SNIP_PDEBC_SYSTEM_5 TOP */
11307                            ///////////////
11308                            // process d //
11309                            ///////////////
11310                            if (add_EM_S) {
11311                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
11312                                if (d.actsExpanded()) {
11313                                    for (index_t k=0; k<numEq; k++) {
11314                                        for (index_t m=0; m<numComp; m++) {
11315                                            const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
11316                                            const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
11317                                            const double d_2 = d_p[INDEX3(k, m, 2, numEq, numComp)];
11318                                            const double d_3 = d_p[INDEX3(k, m, 3, numEq, numComp)];
11319                                            const double tmp2_0 = d_1 + d_3;
11320                                            const double tmp0_0 = d_0 + d_1 + d_2 + d_3;
11321                                            const double tmp1_0 = d_0 + d_2;
11322                                            const double tmp3_0 = d_0 + d_1;
11323                                            const double tmp6_0 = d_0 + d_3;
11324                                            const double tmp4_0 = d_2 + d_3;
11325                                            const double tmp5_0 = d_1 + d_2;
11326                                            const double tmp1_1 = tmp1_0*w25;
11327                                            const double tmp11_1 = d_0*w27;
11328                                            const double tmp9_1 = tmp5_0*w26;
11329                                            const double tmp12_1 = tmp2_0*w25;
11330                                            const double tmp3_1 = tmp3_0*w25;
11331                                            const double tmp6_1 = tmp3_0*w24;
11332                                            const double tmp2_1 = tmp2_0*w24;
11333                                            const double tmp10_1 = d_3*w28;
11334                                            const double tmp16_1 = tmp6_0*w26;
11335                                            const double tmp17_1 = d_1*w28;
11336                                            const double tmp7_1 = d_0*w28;
11337                                            const double tmp14_1 = d_2*w28;
11338                                            const double tmp18_1 = d_2*w27;
11339                                            const double tmp15_1 = d_1*w27;
11340                                            const double tmp0_1 = tmp0_0*w26;
11341                                            const double tmp4_1 = tmp4_0*w24;
11342                                            const double tmp5_1 = tmp4_0*w25;
11343                                            const double tmp8_1 = d_3*w27;
11344                                            const double tmp13_1 = tmp1_0*w24;
11345                                            EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0_1;
11346                                            EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp1_1 + tmp2_1;
11347                                            EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1;
11348                                            EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1;
11349                                            EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp5_1 + tmp6_1;
11350                                            EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp5_1 + tmp6_1;
11351                                            EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp7_1 + tmp8_1 + tmp9_1;
11352                                            EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp9_1;
11353                                            EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp12_1 + tmp13_1;
11354                                            EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp14_1 + tmp15_1 + tmp16_1;
11355                                            EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp3_1 + tmp4_1;
11356                                            EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp16_1 + tmp17_1 + tmp18_1;
11357                                            EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp12_1 + tmp13_1;
11358                                            EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0_1;
11359                                            EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp1_1 + tmp2_1;
11360                                            EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0_1;
11361                                        }
11362                                     }
11363                                } else { /* constant data */
11364                                    for (index_t k=0; k<numEq; k++) {
11365                                        for (index_t m=0; m<numComp; m++) {
11366                                            const double d_0 = d_p[INDEX2(k, m, numEq)];
11367                                            const double tmp2_1 = d_0*w31;
11368                                            const double tmp0_1 = d_0*w30;
11369                                            const double tmp1_1 = d_0*w29;
11370                                            EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0_1;
11371                                            EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp1_1;
11372                                            EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp1_1;
11373                                            EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1;
11374                                            EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp1_1;
11375                                            EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp1_1;
11376                                            EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp2_1;
11377                                            EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp2_1;
11378                                            EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp1_1;
11379                                            EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp2_1;
11380                                            EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp1_1;
11381                                            EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp2_1;
11382                                            EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp1_1;
11383                                            EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0_1;
11384                                            EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp1_1;
11385                                            EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0_1;
11386                                        }
11387                                    }
11388                                }
11389                            }
11390                            ///////////////
11391                            // process y //
11392                            ///////////////
11393                            if (add_EM_F) {
11394                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
11395                                if (y.actsExpanded()) {
11396                                    for (index_t k=0; k<numEq; k++) {
11397                                        const double y_0 = y_p[INDEX2(k, 0, numEq)];
11398                                        const double y_1 = y_p[INDEX2(k, 1, numEq)];
11399                                        const double y_2 = y_p[INDEX2(k, 2, numEq)];
11400                                        const double y_3 = y_p[INDEX2(k, 3, numEq)];
11401                                        const double tmp0_0 = y_1 + y_2;
11402                                        const double tmp1_0 = y_0 + y_3;
11403                                        const double tmp7_1 = w34*y_1;
11404                                        const double tmp3_1 = w32*y_1;
11405                                        const double tmp8_1 = w32*y_3;
11406                                        const double tmp0_1 = w32*y_0;
11407                                        const double tmp2_1 = w34*y_3;
11408                                        const double tmp9_1 = w34*y_0;
11409                                        const double tmp6_1 = w32*y_2;
11410                                        const double tmp5_1 = w34*y_2;
11411                                        const double tmp1_1 = tmp0_0*w33;
11412                                        const double tmp4_1 = tmp1_0*w33;
11413                                        EM_F[INDEX2(k,4,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
11414                                        EM_F[INDEX2(k,5,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
11415                                        EM_F[INDEX2(k,6,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
11416                                        EM_F[INDEX2(k,7,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;
11417                                    }
11418                                } else { /* constant data */
11419                                    for (index_t k=0; k<numEq; k++) {
11420                                        const double y_0 = y_p[k];
11421                                        const double tmp0_1 = w35*y_0;
11422                                        EM_F[INDEX2(k,4,numEq)]+=tmp0_1;
11423                                        EM_F[INDEX2(k,5,numEq)]+=tmp0_1;
11424                                        EM_F[INDEX2(k,6,numEq)]+=tmp0_1;
11425                                        EM_F[INDEX2(k,7,numEq)]+=tmp0_1;
11426                                    }
11427                                }
11428                            }
11429                            /* GENERATOR SNIP_PDEBC_SYSTEM_5 BOTTOM */
11430                            const index_t firstNode=m_N0*m_N1*(m_N2-2)+m_N0*k1+k0;
11431                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11432                                    add_EM_F, firstNode, numEq, numComp);
11433                        } // k0 loop
11434                    } // k1 loop
11435                } // colouring
11436            } // face 5
11437        } // end of parallel region
11438  }  }
11439    
11440  //protected  //protected
11441  void Brick::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,  void Brick::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
11442        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
11443  {  {
11444        const double h0 = m_l0/m_gNE0;
11445        const double h1 = m_l1/m_gNE1;
11446        const double h2 = m_l2/m_gNE2;
11447        dim_t numEq, numComp;
11448        if (!mat)
11449            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
11450        else {
11451            numEq=mat->logical_row_block_size;
11452            numComp=mat->logical_col_block_size;
11453        }
11454        /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_PRE TOP */
11455        const double w0 = 0.0625*h1*h2;
11456        const double w1 = 0.25*h1*h2;
11457        const double w2 = 0.0625*h0*h2;
11458        const double w3 = 0.25*h0*h2;
11459        const double w4 = 0.0625*h0*h1;
11460        const double w5 = 0.25*h0*h1;
11461        /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_PRE BOTTOM */
11462        const bool add_EM_S=!d.isEmpty();
11463        const bool add_EM_F=!y.isEmpty();
11464        rhs.requireWrite();
11465    #pragma omp parallel
11466        {
11467            if (m_faceOffset[0] > -1) {
11468                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11469    #pragma omp for nowait
11470                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11471                        for (index_t k1=0; k1<m_NE1; ++k1) {
11472                            vector<double> EM_S(8*8*numEq*numComp, 0);
11473                            vector<double> EM_F(8*numEq, 0);
11474                            const index_t e = INDEX2(k1,k2,m_NE1);
11475                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_0 TOP */
11476                            ///////////////
11477                            // process d //
11478                            ///////////////
11479                            if (add_EM_S) {
11480                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
11481                                for (index_t k=0; k<numEq; k++) {
11482                                    for (index_t m=0; m<numComp; m++) {
11483                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
11484                                        const double tmp0_1 = d_0*w0;
11485                                        EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0_1;
11486                                        EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1;
11487                                        EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0_1;
11488                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0_1;
11489                                        EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0_1;
11490                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0_1;
11491                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0_1;
11492                                        EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0_1;
11493                                        EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0_1;
11494                                        EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp0_1;
11495                                        EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp0_1;
11496                                        EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1;
11497                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1;
11498                                        EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp0_1;
11499                                        EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1;
11500                                        EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp0_1;
11501                                    }
11502                                }
11503                            }
11504                            ///////////////
11505                            // process y //
11506                            ///////////////
11507                            if (add_EM_F) {
11508                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
11509                                for (index_t k=0; k<numEq; k++) {
11510                                    const double y_0 = y_p[k];
11511                                    const double tmp0_1 = w1*y_0;
11512                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
11513                                    EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
11514                                    EM_F[INDEX2(k,4,numEq)]+=tmp0_1;
11515                                    EM_F[INDEX2(k,6,numEq)]+=tmp0_1;
11516                                }
11517                            }
11518                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_0 BOTTOM */
11519                            const index_t firstNode=m_N0*m_N1*k2+m_N0*k1;
11520                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11521                                    add_EM_F, firstNode, numEq, numComp);
11522                        } // k1 loop
11523                    } // k2 loop
11524                } // colouring
11525            } // face 0
11526    
11527            if (m_faceOffset[1] > -1) {
11528                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11529    #pragma omp for nowait
11530                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11531                        for (index_t k1=0; k1<m_NE1; ++k1) {
11532                            vector<double> EM_S(8*8*numEq*numComp, 0);
11533                            vector<double> EM_F(8*numEq, 0);
11534                            const index_t e = m_faceOffset[1]+INDEX2(k1,k2,m_NE1);
11535                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_1 TOP */
11536                            ///////////////
11537                            // process d //
11538                            ///////////////
11539                            if (add_EM_S) {
11540                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
11541                                for (index_t k=0; k<numEq; k++) {
11542                                    for (index_t m=0; m<numComp; m++) {
11543                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
11544                                        const double tmp0_1 = d_0*w0;
11545                                        EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1;
11546                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1;
11547                                        EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0_1;
11548                                        EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1;
11549                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0_1;
11550                                        EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0_1;
11551                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0_1;
11552                                        EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0_1;
11553                                        EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp0_1;
11554                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0_1;
11555                                        EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp0_1;
11556                                        EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0_1;
11557                                        EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0_1;
11558                                        EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp0_1;
11559                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1;
11560                                        EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp0_1;
11561                                    }
11562                                }
11563                            }
11564                            ///////////////
11565                            // process y //
11566                            ///////////////
11567                            if (add_EM_F) {
11568                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
11569                                for (index_t k=0; k<numEq; k++) {
11570                                    const double y_0 = y_p[k];
11571                                    const double tmp0_1 = w1*y_0;
11572                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
11573                                    EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
11574                                    EM_F[INDEX2(k,5,numEq)]+=tmp0_1;
11575                                    EM_F[INDEX2(k,7,numEq)]+=tmp0_1;
11576                                }
11577                            }
11578                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_1 BOTTOM */
11579                            const index_t firstNode=m_N0*m_N1*k2+m_N0*(k1+1)-2;
11580                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11581                                    add_EM_F, firstNode, numEq, numComp);
11582                        } // k1 loop
11583                    } // k2 loop
11584                } // colouring
11585            } // face 1
11586    
11587            if (m_faceOffset[2] > -1) {
11588                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11589    #pragma omp for nowait
11590                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11591                        for (index_t k0=0; k0<m_NE0; ++k0) {
11592                            vector<double> EM_S(8*8*numEq*numComp, 0);
11593                            vector<double> EM_F(8*numEq, 0);
11594                            const index_t e = m_faceOffset[2]+INDEX2(k0,k2,m_NE0);
11595                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_2 TOP */
11596                            ///////////////
11597                            // process d //
11598                            ///////////////
11599                            if (add_EM_S) {
11600                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
11601                                for (index_t k=0; k<numEq; k++) {
11602                                    for (index_t m=0; m<numComp; m++) {
11603                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
11604                                        const double tmp0_1 = d_0*w2;
11605                                        EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0_1;
11606                                        EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0_1;
11607                                        EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0_1;
11608                                        EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0_1;
11609                                        EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0_1;
11610                                        EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp0_1;
11611                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0_1;
11612                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0_1;
11613                                        EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp0_1;
11614                                        EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0_1;
11615                                        EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0_1;
11616                                        EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0_1;
11617                                        EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp0_1;
11618                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0_1;
11619                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0_1;
11620                                        EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp0_1;
11621                                    }
11622                                }
11623                            }
11624                            ///////////////
11625                            // process y //
11626                            ///////////////
11627                            if (add_EM_F) {
11628                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
11629                                for (index_t k=0; k<numEq; k++) {
11630                                    const double y_0 = y_p[k];
11631                                    const double tmp0_1 = w3*y_0;
11632                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
11633                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
11634                                    EM_F[INDEX2(k,4,numEq)]+=tmp0_1;
11635                                    EM_F[INDEX2(k,5,numEq)]+=tmp0_1;
11636                                }
11637                            }
11638                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_2 BOTTOM */
11639                            const index_t firstNode=m_N0*m_N1*k2+k0;
11640                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11641                                    add_EM_F, firstNode, numEq, numComp);
11642                        } // k0 loop
11643                    } // k2 loop
11644                } // colouring
11645            } // face 2
11646    
11647            if (m_faceOffset[3] > -1) {
11648                for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11649    #pragma omp for nowait
11650                    for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11651                        for (index_t k0=0; k0<m_NE0; ++k0) {
11652                            vector<double> EM_S(8*8*numEq*numComp, 0);
11653                            vector<double> EM_F(8*numEq, 0);
11654                            const index_t e = m_faceOffset[3]+INDEX2(k0,k2,m_NE0);
11655                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_3 TOP */
11656                            ///////////////
11657                            // process d //
11658                            ///////////////
11659                            if (add_EM_S) {
11660                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
11661                                for (index_t k=0; k<numEq; k++) {
11662                                    for (index_t m=0; m<numComp; m++) {
11663                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
11664                                        const double tmp0_1 = d_0*w2;
11665                                        EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1;
11666                                        EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1;
11667                                        EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1;
11668                                        EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0_1;
11669                                        EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0_1;
11670                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0_1;
11671                                        EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0_1;
11672                                        EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0_1;
11673                                        EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0_1;
11674                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0_1;
11675                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0_1;
11676                                        EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0_1;
11677                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1;
11678                                        EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0_1;
11679                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1;
11680                                        EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1;
11681                                    }
11682                                }
11683                            }
11684                            ///////////////
11685                            // process y //
11686                            ///////////////
11687                            if (add_EM_F) {
11688                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
11689                                for (index_t k=0; k<numEq; k++) {
11690                                    const double y_0 = y_p[k];
11691                                    const double tmp0_1 = w3*y_0;
11692                                    EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
11693                                    EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
11694                                    EM_F[INDEX2(k,6,numEq)]+=tmp0_1;
11695                                    EM_F[INDEX2(k,7,numEq)]+=tmp0_1;
11696                                }
11697                            }
11698                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_3 BOTTOM */
11699                            const index_t firstNode=m_N0*m_N1*k2+m_N0*(m_N1-2)+k0;
11700                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11701                                    add_EM_F, firstNode, numEq, numComp);
11702                        } // k0 loop
11703                    } // k2 loop
11704                } // colouring
11705            } // face 3
11706    
11707            if (m_faceOffset[4] > -1) {
11708                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
11709    #pragma omp for nowait
11710                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
11711                        for (index_t k0=0; k0<m_NE0; ++k0) {
11712                            vector<double> EM_S(8*8*numEq*numComp, 0);
11713                            vector<double> EM_F(8*numEq, 0);
11714                            const index_t e = m_faceOffset[4]+INDEX2(k0,k1,m_NE0);
11715                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_4 TOP */
11716                            ///////////////
11717                            // process d //
11718                            ///////////////
11719                            if (add_EM_S) {
11720                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
11721                                for (index_t k=0; k<numEq; k++) {
11722                                    for (index_t m=0; m<numComp; m++) {
11723                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
11724                                        const double tmp0_1 = d_0*w4;
11725                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1;
11726                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp0_1;
11727                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp0_1;
11728                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0_1;
11729                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0_1;
11730                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0_1;
11731                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0_1;
11732                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0_1;
11733                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1;
11734                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0_1;
11735                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1;
11736                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0_1;
11737                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0_1;
11738                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0_1;
11739                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1;
11740                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1;
11741                                    }
11742                                }
11743                            }
11744                            ///////////////
11745                            // process y //
11746                            ///////////////
11747                            if (add_EM_F) {
11748                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
11749                                for (index_t k=0; k<numEq; k++) {
11750                                    const double y_0 = y_p[k];
11751                                    const double tmp0_1 = w5*y_0;
11752                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
11753                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
11754                                    EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
11755                                    EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
11756                                }
11757                            }
11758                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_4 BOTTOM */
11759                            const index_t firstNode=m_N0*k1+k0;
11760                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11761                                    add_EM_F, firstNode, numEq, numComp);
11762                        } // k0 loop
11763                    } // k1 loop
11764                } // colouring
11765            } // face 4
11766    
11767            if (m_faceOffset[5] > -1) {
11768                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
11769    #pragma omp for nowait
11770                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
11771                        for (index_t k0=0; k0<m_NE0; ++k0) {
11772                            vector<double> EM_S(8*8*numEq*numComp, 0);
11773                            vector<double> EM_F(8*numEq, 0);
11774                            const index_t e = m_faceOffset[5]+INDEX2(k0,k1,m_NE0);
11775                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_5 TOP */
11776                            ///////////////
11777                            // process d //
11778                            ///////////////
11779                            if (add_EM_S) {
11780                                const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
11781                                for (index_t k=0; k<numEq; k++) {
11782                                    for (index_t m=0; m<numComp; m++) {
11783                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
11784                                        const double tmp0_1 = d_0*w4;
11785                                        EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0_1;
11786                                        EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0_1;
11787                                        EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0_1;
11788                                        EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1;
11789                                        EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0_1;
11790                                        EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0_1;
11791                                        EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0_1;
11792                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0_1;
11793                                        EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0_1;
11794                                        EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0_1;
11795                                        EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0_1;
11796                                        EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0_1;
11797                                        EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0_1;
11798                                        EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0_1;
11799                                        EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1;
11800                                        EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0_1;
11801                                    }
11802                                }
11803                            }
11804                            ///////////////
11805                            // process y //
11806                            ///////////////
11807                            if (add_EM_F) {
11808                                const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
11809                                for (index_t k=0; k<numEq; k++) {
11810                                    const double y_0 = y_p[k];
11811                                    const double tmp0_1 = w5*y_0;
11812                                    EM_F[INDEX2(k,4,numEq)]+=tmp0_1;
11813                                    EM_F[INDEX2(k,5,numEq)]+=tmp0_1;
11814                                    EM_F[INDEX2(k,6,numEq)]+=tmp0_1;
11815                                    EM_F[INDEX2(k,7,numEq)]+=tmp0_1;
11816                                }
11817                            }
11818                            /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_5 BOTTOM */
11819                            const index_t firstNode=m_N0*m_N1*(m_N2-2)+m_N0*k1+k0;
11820                            addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S,
11821                                    add_EM_F, firstNode, numEq, numComp);
11822                        } // k0 loop
11823                    } // k1 loop
11824                } // colouring
11825            } // face 5
11826        } // end of parallel region
11827  }  }
11828    
11829    

Legend:
Removed from v.3769  
changed lines
  Added in v.3777

  ViewVC Help
Powered by ViewVC 1.1.26