/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

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

revision 3777 by caltinay, Thu Jan 19 06:17:38 2012 UTC revision 3791 by caltinay, Wed Feb 1 05:10:22 2012 UTC
# Line 13  Line 13 
13    
14  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
15  extern "C" {  extern "C" {
16  #include "paso/SystemMatrix.h"  #include <paso/SystemMatrix.h>
17  }  }
18    
19  #if USE_SILO  #if USE_SILO
# Line 29  using namespace std; Line 29  using namespace std;
29    
30  namespace ripley {  namespace ripley {
31    
32  Rectangle::Rectangle(int n0, int n1, double l0, double l1, int d0, int d1) :  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
33                         double y1, int d0, int d1) :
34      RipleyDomain(2),      RipleyDomain(2),
35      m_gNE0(n0),      m_gNE0(n0),
36      m_gNE1(n1),      m_gNE1(n1),
37      m_l0(l0),      m_x0(x0),
38      m_l1(l1),      m_y0(y0),
39        m_l0(x1-x0),
40        m_l1(y1-y0),
41      m_NX(d0),      m_NX(d0),
42      m_NY(d1)      m_NY(d1)
43  {  {
# Line 80  Rectangle::Rectangle(int n0, int n1, dou Line 83  Rectangle::Rectangle(int n0, int n1, dou
83    
84  Rectangle::~Rectangle()  Rectangle::~Rectangle()
85  {  {
86        Paso_SystemMatrixPattern_free(m_pattern);
87        Paso_Connector_free(m_connector);
88  }  }
89    
90  string Rectangle::getDescription() const  string Rectangle::getDescription() const
# Line 93  bool Rectangle::operator==(const Abstrac Line 98  bool Rectangle::operator==(const Abstrac
98      if (o) {      if (o) {
99          return (RipleyDomain::operator==(other) &&          return (RipleyDomain::operator==(other) &&
100                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1
101                    && m_x0==o->m_x0 && m_y0==o->m_y0
102                  && m_l0==o->m_l0 && m_l1==o->m_l1                  && m_l0==o->m_l0 && m_l1==o->m_l1
103                  && m_NX==o->m_NX && m_NY==o->m_NY);                  && m_NX==o->m_NX && m_NY==o->m_NY);
104      }      }
# Line 240  void Rectangle::dump(const string& fileN Line 246  void Rectangle::dump(const string& fileN
246      }      }
247    
248  #else // USE_SILO  #else // USE_SILO
249      throw RipleyException("dump(): no Silo support");      throw RipleyException("dump: no Silo support");
250  #endif  #endif
251  }  }
252    
# Line 264  const int* Rectangle::borrowSampleRefere Line 270  const int* Rectangle::borrowSampleRefere
270      }      }
271    
272      stringstream msg;      stringstream msg;
273      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
274      throw RipleyException(msg.str());      throw RipleyException(msg.str());
275  }  }
276    
# Line 314  bool Rectangle::ownSample(int fsType, in Line 319  bool Rectangle::ownSample(int fsType, in
319      }      }
320    
321      stringstream msg;      stringstream msg;
322      msg << "ownSample() not implemented for "      msg << "ownSample: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
323      throw RipleyException(msg.str());      throw RipleyException(msg.str());
324  }  }
325    
# Line 416  void Rectangle::setToNormal(escript::Dat Line 420  void Rectangle::setToNormal(escript::Dat
420    
421      } else {      } else {
422          stringstream msg;          stringstream msg;
423          msg << "setToNormal() not implemented for "          msg << "setToNormal: invalid function space type "
424              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());              << out.getFunctionSpace().getTypeCode();
425          throw RipleyException(msg.str());          throw RipleyException(msg.str());
426      }      }
427  }  }
# Line 479  void Rectangle::setToSize(escript::Data& Line 483  void Rectangle::setToSize(escript::Data&
483    
484      } else {      } else {
485          stringstream msg;          stringstream msg;
486          msg << "setToSize() not implemented for "          msg << "setToSize: invalid function space type "
487              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());              << out.getFunctionSpace().getTypeCode();
488          throw RipleyException(msg.str());          throw RipleyException(msg.str());
489      }      }
490  }  }
# Line 557  IndexVector Rectangle::getNumSubdivision Line 561  IndexVector Rectangle::getNumSubdivision
561  pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
562  {  {
563      if (dim==0) {      if (dim==0) {
564          return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);          return pair<double,double>(m_x0+(m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
565      } else if (dim==1) {      } else if (dim==1) {
566          return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);          return pair<double,double>(m_y0+(m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
567      }      }
568      throw RipleyException("getFirstCoordAndSpacing(): invalid argument");      throw RipleyException("getFirstCoordAndSpacing: invalid argument");
569  }  }
570    
571  //protected  //protected
# Line 2725  void Rectangle::assemblePDESystem(Paso_S Line 2729  void Rectangle::assemblePDESystem(Paso_S
2729                                      const double tmp78_1 = B_0_0*w34;                                      const double tmp78_1 = B_0_0*w34;
2730                                      const double tmp12_1 = tmp0_0*w36;                                      const double tmp12_1 = tmp0_0*w36;
2731                                      const double tmp75_1 = B_1_0*w32;                                      const double tmp75_1 = B_1_0*w32;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
2732                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
2733                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2734                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2735                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
2736                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2737                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
2738                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
2739                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
2740                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
2741                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
2742                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
2743                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2744                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
2745                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
2746                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2747                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2748                                  }                                  }
2749                              }                              }
2750                          } else { // constant data                          } else { // constant data
# Line 2756  void Rectangle::assemblePDESystem(Paso_S Line 2760  void Rectangle::assemblePDESystem(Paso_S
2760                                      const double tmp2_1 = B_0*w46;                                      const double tmp2_1 = B_0*w46;
2761                                      const double tmp7_1 = B_0*w51;                                      const double tmp7_1 = B_0*w51;
2762                                      const double tmp3_1 = B_0*w47;                                      const double tmp3_1 = B_0*w47;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
2763                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;
2764                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
2765                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;
2766                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
2767                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
2768                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;
2769                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;
2770                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;
2771                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;
2772                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;  
2773                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;
2774                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
2775                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;
2776                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;
2777                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;
2778                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;
2779                                  }                                  }
2780                              }                              }
2781                          }                          }
# Line 3349  void Rectangle::assemblePDEBoundarySingl Line 3353  void Rectangle::assemblePDEBoundarySingl
3353  {  {
3354      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
3355      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
     /* GENERATOR SNIP_PDEBC_SINGLE_PRE TOP */  
3356      const double w0 = 0.31100423396407310779*h1;      const double w0 = 0.31100423396407310779*h1;
3357      const double w1 = 0.022329099369260225539*h1;      const double w1 = 0.022329099369260225539*h1;
3358      const double w10 = 0.022329099369260225539*h0;      const double w10 = 0.022329099369260225539*h0;
# Line 3366  void Rectangle::assemblePDEBoundarySingl Line 3369  void Rectangle::assemblePDEBoundarySingl
3369      const double w7 = 0.5*h1;      const double w7 = 0.5*h1;
3370      const double w8 = 0.083333333333333333333*h0;      const double w8 = 0.083333333333333333333*h0;
3371      const double w9 = 0.31100423396407310779*h0;      const double w9 = 0.31100423396407310779*h0;
     /* GENERATOR SNIP_PDEBC_SINGLE_PRE BOTTOM */  
3372      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3373      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3374      rhs.requireWrite();      rhs.requireWrite();
# Line 3379  void Rectangle::assemblePDEBoundarySingl Line 3381  void Rectangle::assemblePDEBoundarySingl
3381                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3382                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3383                      const index_t e = k1;                      const index_t e = k1;
                     /* GENERATOR SNIP_PDEBC_SINGLE_0 TOP */  
3384                      ///////////////                      ///////////////
3385                      // process d //                      // process d //
3386                      ///////////////                      ///////////////
# Line 3395  void Rectangle::assemblePDEBoundarySingl Line 3396  void Rectangle::assemblePDEBoundarySingl
3396                              const double tmp3_1 = d_1*w0;                              const double tmp3_1 = d_1*w0;
3397                              const double tmp2_1 = tmp0_0*w2;                              const double tmp2_1 = tmp0_0*w2;
3398                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
                             EM_S[INDEX2(0,2,4)]+=tmp2_1;  
3399                              EM_S[INDEX2(2,0,4)]+=tmp2_1;                              EM_S[INDEX2(2,0,4)]+=tmp2_1;
3400                                EM_S[INDEX2(0,2,4)]+=tmp2_1;
3401                              EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;
3402                          } else { /* constant data */                          } else { /* constant data */
3403                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3404                              const double tmp1_1 = d_0*w4;                              const double tmp1_1 = d_0*w4;
3405                              const double tmp0_1 = d_0*w3;                              const double tmp0_1 = d_0*w3;
3406                              EM_S[INDEX2(0,0,4)]+=tmp0_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1;
                             EM_S[INDEX2(0,2,4)]+=tmp1_1;  
3407                              EM_S[INDEX2(2,0,4)]+=tmp1_1;                              EM_S[INDEX2(2,0,4)]+=tmp1_1;
3408                                EM_S[INDEX2(0,2,4)]+=tmp1_1;
3409                              EM_S[INDEX2(2,2,4)]+=tmp0_1;                              EM_S[INDEX2(2,2,4)]+=tmp0_1;
3410                          }                          }
3411                      }                      }
# Line 3429  void Rectangle::assemblePDEBoundarySingl Line 3430  void Rectangle::assemblePDEBoundarySingl
3430                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
3431                          }                          }
3432                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_0 BOTTOM */  
3433                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_N0*k1;
3434                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3435                  }                  }
# Line 3443  void Rectangle::assemblePDEBoundarySingl Line 3443  void Rectangle::assemblePDEBoundarySingl
3443                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3444                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3445                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
                     /* GENERATOR SNIP_PDEBC_SINGLE_1 TOP */  
3446                      ///////////////                      ///////////////
3447                      // process d //                      // process d //
3448                      ///////////////                      ///////////////
# Line 3458  void Rectangle::assemblePDEBoundarySingl Line 3457  void Rectangle::assemblePDEBoundarySingl
3457                              const double tmp3_1 = d_0*w0;                              const double tmp3_1 = d_0*w0;
3458                              const double tmp0_1 = d_1*w0;                              const double tmp0_1 = d_1*w0;
3459                              const double tmp2_1 = tmp0_0*w2;                              const double tmp2_1 = tmp0_0*w2;
3460                              EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;
3461                              EM_S[INDEX2(3,1,4)]+=tmp2_1;                              EM_S[INDEX2(3,1,4)]+=tmp2_1;
3462                              EM_S[INDEX2(1,3,4)]+=tmp2_1;                              EM_S[INDEX2(1,3,4)]+=tmp2_1;
3463                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp1_1;
3464                          } else { /* constant data */                          } else { /* constant data */
3465                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3466                              const double tmp1_1 = d_0*w4;                              const double tmp1_1 = d_0*w4;
3467                              const double tmp0_1 = d_0*w3;                              const double tmp0_1 = d_0*w3;
3468                              EM_S[INDEX2(3,3,4)]+=tmp0_1;                              EM_S[INDEX2(1,1,4)]+=tmp0_1;
3469                              EM_S[INDEX2(3,1,4)]+=tmp1_1;                              EM_S[INDEX2(3,1,4)]+=tmp1_1;
3470                              EM_S[INDEX2(1,3,4)]+=tmp1_1;                              EM_S[INDEX2(1,3,4)]+=tmp1_1;
3471                              EM_S[INDEX2(1,1,4)]+=tmp0_1;                              EM_S[INDEX2(3,3,4)]+=tmp0_1;
3472                          }                          }
3473                      }                      }
3474                      ///////////////                      ///////////////
# Line 3493  void Rectangle::assemblePDEBoundarySingl Line 3492  void Rectangle::assemblePDEBoundarySingl
3492                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
3493                          }                          }
3494                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_1 BOTTOM */  
3495                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_N0*(k1+1)-2;
3496                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3497                  }                  }
# Line 3507  void Rectangle::assemblePDEBoundarySingl Line 3505  void Rectangle::assemblePDEBoundarySingl
3505                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3506                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3507                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
                     /* GENERATOR SNIP_PDEBC_SINGLE_2 TOP */  
3508                      ///////////////                      ///////////////
3509                      // process d //                      // process d //
3510                      ///////////////                      ///////////////
# Line 3522  void Rectangle::assemblePDEBoundarySingl Line 3519  void Rectangle::assemblePDEBoundarySingl
3519                              const double tmp0_1 = tmp0_0*w8;                              const double tmp0_1 = tmp0_0*w8;
3520                              const double tmp1_1 = d_1*w10;                              const double tmp1_1 = d_1*w10;
3521                              const double tmp3_1 = d_0*w10;                              const double tmp3_1 = d_0*w10;
                             EM_S[INDEX2(0,1,4)]+=tmp0_1;  
3522                              EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;                              EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;
3523                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
3524                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
3525                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;
3526                          } else { /* constant data */                          } else { /* constant data */
3527                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3528                              const double tmp0_1 = d_0*w11;                              const double tmp0_1 = d_0*w11;
3529                              const double tmp1_1 = d_0*w12;                              const double tmp1_1 = d_0*w12;
                             EM_S[INDEX2(0,1,4)]+=tmp0_1;  
3530                              EM_S[INDEX2(0,0,4)]+=tmp1_1;                              EM_S[INDEX2(0,0,4)]+=tmp1_1;
3531                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
3532                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
3533                              EM_S[INDEX2(1,1,4)]+=tmp1_1;                              EM_S[INDEX2(1,1,4)]+=tmp1_1;
3534                          }                          }
3535                      }                      }
# Line 3557  void Rectangle::assemblePDEBoundarySingl Line 3554  void Rectangle::assemblePDEBoundarySingl
3554                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
3555                          }                          }
3556                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_2 BOTTOM */  
3557                      const index_t firstNode=k0;                      const index_t firstNode=k0;
3558                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3559                  }                  }
# Line 3571  void Rectangle::assemblePDEBoundarySingl Line 3567  void Rectangle::assemblePDEBoundarySingl
3567                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
3568                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3569                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
                     /* GENERATOR SNIP_PDEBC_SINGLE_3 TOP */  
3570                      ///////////////                      ///////////////
3571                      // process d //                      // process d //
3572                      ///////////////                      ///////////////
# Line 3586  void Rectangle::assemblePDEBoundarySingl Line 3581  void Rectangle::assemblePDEBoundarySingl
3581                              const double tmp0_1 = tmp0_0*w8;                              const double tmp0_1 = tmp0_0*w8;
3582                              const double tmp3_1 = d_1*w10;                              const double tmp3_1 = d_1*w10;
3583                              const double tmp1_1 = d_0*w10;                              const double tmp1_1 = d_0*w10;
3584                                EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;
3585                              EM_S[INDEX2(3,2,4)]+=tmp0_1;                              EM_S[INDEX2(3,2,4)]+=tmp0_1;
                             EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;  
3586                              EM_S[INDEX2(2,3,4)]+=tmp0_1;                              EM_S[INDEX2(2,3,4)]+=tmp0_1;
3587                              EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
3588                          } else { /* constant data */                          } else { /* constant data */
3589                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3590                              const double tmp0_1 = d_0*w11;                              const double tmp0_1 = d_0*w11;
3591                              const double tmp1_1 = d_0*w12;                              const double tmp1_1 = d_0*w12;
3592                                EM_S[INDEX2(2,2,4)]+=tmp1_1;
3593                              EM_S[INDEX2(3,2,4)]+=tmp0_1;                              EM_S[INDEX2(3,2,4)]+=tmp0_1;
                             EM_S[INDEX2(3,3,4)]+=tmp1_1;  
3594                              EM_S[INDEX2(2,3,4)]+=tmp0_1;                              EM_S[INDEX2(2,3,4)]+=tmp0_1;
3595                              EM_S[INDEX2(2,2,4)]+=tmp1_1;                              EM_S[INDEX2(3,3,4)]+=tmp1_1;
3596                          }                          }
3597                      }                      }
3598                      ///////////////                      ///////////////
# Line 3621  void Rectangle::assemblePDEBoundarySingl Line 3616  void Rectangle::assemblePDEBoundarySingl
3616                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
3617                          }                          }
3618                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_3 BOTTOM */  
3619                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_N0*(m_N1-2)+k0;
3620                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3621                  }                  }
# Line 3636  void Rectangle::assemblePDEBoundarySingl Line 3630  void Rectangle::assemblePDEBoundarySingl
3630  {  {
3631      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
3632      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE TOP */  
3633      const double w0 = 0.25*h1;      const double w0 = 0.25*h1;
3634      const double w1 = 0.5*h1;      const double w1 = 0.5*h1;
3635      const double w2 = 0.25*h0;      const double w2 = 0.25*h0;
3636      const double w3 = 0.5*h0;      const double w3 = 0.5*h0;
     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE BOTTOM */  
3637      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3638      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3639      rhs.requireWrite();      rhs.requireWrite();
# Line 3654  void Rectangle::assemblePDEBoundarySingl Line 3646  void Rectangle::assemblePDEBoundarySingl
3646                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3647                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3648                      const index_t e = k1;                      const index_t e = k1;
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 TOP */  
3649                      ///////////////                      ///////////////
3650                      // process d //                      // process d //
3651                      ///////////////                      ///////////////
# Line 3663  void Rectangle::assemblePDEBoundarySingl Line 3654  void Rectangle::assemblePDEBoundarySingl
3654                          const double d_0 = d_p[0];                          const double d_0 = d_p[0];
3655                          const double tmp0_1 = d_0*w0;                          const double tmp0_1 = d_0*w0;
3656                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1;
                         EM_S[INDEX2(0,2,4)]+=tmp0_1;  
3657                          EM_S[INDEX2(2,0,4)]+=tmp0_1;                          EM_S[INDEX2(2,0,4)]+=tmp0_1;
3658                            EM_S[INDEX2(0,2,4)]+=tmp0_1;
3659                          EM_S[INDEX2(2,2,4)]+=tmp0_1;                          EM_S[INDEX2(2,2,4)]+=tmp0_1;
3660                      }                      }
3661                      ///////////////                      ///////////////
# Line 3672  void Rectangle::assemblePDEBoundarySingl Line 3663  void Rectangle::assemblePDEBoundarySingl
3663                      ///////////////                      ///////////////
3664                      if (add_EM_F) {                      if (add_EM_F) {
3665                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3666                          const double y_0 = y_p[0];                          const double tmp0_1 = w1*y_p[0];
                         const double tmp0_1 = w1*y_0;  
3667                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0_1;
3668                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
3669                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 BOTTOM */  
3670                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_N0*k1;
3671                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3672                  }                  }
# Line 3691  void Rectangle::assemblePDEBoundarySingl Line 3680  void Rectangle::assemblePDEBoundarySingl
3680                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3681                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3682                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 TOP */  
3683                      ///////////////                      ///////////////
3684                      // process d //                      // process d //
3685                      ///////////////                      ///////////////
# Line 3699  void Rectangle::assemblePDEBoundarySingl Line 3687  void Rectangle::assemblePDEBoundarySingl
3687                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3688                          const double d_0 = d_p[0];                          const double d_0 = d_p[0];
3689                          const double tmp0_1 = d_0*w0;                          const double tmp0_1 = d_0*w0;
3690                          EM_S[INDEX2(3,3,4)]+=tmp0_1;                          EM_S[INDEX2(1,1,4)]+=tmp0_1;
3691                          EM_S[INDEX2(3,1,4)]+=tmp0_1;                          EM_S[INDEX2(3,1,4)]+=tmp0_1;
3692                          EM_S[INDEX2(1,3,4)]+=tmp0_1;                          EM_S[INDEX2(1,3,4)]+=tmp0_1;
3693                          EM_S[INDEX2(1,1,4)]+=tmp0_1;                          EM_S[INDEX2(3,3,4)]+=tmp0_1;
3694                      }                      }
3695                      ///////////////                      ///////////////
3696                      // process y //                      // process y //
3697                      ///////////////                      ///////////////
3698                      if (add_EM_F) {                      if (add_EM_F) {
3699                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3700                          const double y_0 = y_p[0];                          const double tmp0_1 = w1*y_p[0];
                         const double tmp0_1 = w1*y_0;  
3701                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0_1;
3702                          EM_F[3]+=tmp0_1;                          EM_F[3]+=tmp0_1;
3703                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 BOTTOM */  
3704                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_N0*(k1+1)-2;
3705                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3706                  }                  }
# Line 3728  void Rectangle::assemblePDEBoundarySingl Line 3714  void Rectangle::assemblePDEBoundarySingl
3714                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3715                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3716                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 TOP */  
3717                      ///////////////                      ///////////////
3718                      // process d //                      // process d //
3719                      ///////////////                      ///////////////
3720                      if (add_EM_S) {                      if (add_EM_S) {
3721                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3722                          const double d_0 = d_p[0];                          const double tmp0_1 = d_p[0]*w2;
                         const double tmp0_1 = d_0*w2;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1;  
3723                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1;
3724                          EM_S[INDEX2(1,0,4)]+=tmp0_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1;
3725                            EM_S[INDEX2(0,1,4)]+=tmp0_1;
3726                          EM_S[INDEX2(1,1,4)]+=tmp0_1;                          EM_S[INDEX2(1,1,4)]+=tmp0_1;
3727                      }                      }
3728                      ///////////////                      ///////////////
# Line 3746  void Rectangle::assemblePDEBoundarySingl Line 3730  void Rectangle::assemblePDEBoundarySingl
3730                      ///////////////                      ///////////////
3731                      if (add_EM_F) {                      if (add_EM_F) {
3732                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3733                          const double y_0 = y_p[0];                          const double tmp0_1 = w3*y_p[0];
                         const double tmp0_1 = w3*y_0;  
3734                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0_1;
3735                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0_1;
3736                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 BOTTOM */  
3737                      const index_t firstNode=k0;                      const index_t firstNode=k0;
3738                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3739                  }                  }
# Line 3765  void Rectangle::assemblePDEBoundarySingl Line 3747  void Rectangle::assemblePDEBoundarySingl
3747                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3748                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3749                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 TOP */  
3750                      ///////////////                      ///////////////
3751                      // process d //                      // process d //
3752                      ///////////////                      ///////////////
3753                      if (add_EM_S) {                      if (add_EM_S) {
3754                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3755                          const double d_0 = d_p[0];                          const double tmp0_1 = d_p[0]*w2;
3756                          const double tmp0_1 = d_0*w2;                          EM_S[INDEX2(2,2,4)]+=tmp0_1;
3757                          EM_S[INDEX2(3,2,4)]+=tmp0_1;                          EM_S[INDEX2(3,2,4)]+=tmp0_1;
                         EM_S[INDEX2(3,3,4)]+=tmp0_1;  
3758                          EM_S[INDEX2(2,3,4)]+=tmp0_1;                          EM_S[INDEX2(2,3,4)]+=tmp0_1;
3759                          EM_S[INDEX2(2,2,4)]+=tmp0_1;                          EM_S[INDEX2(3,3,4)]+=tmp0_1;
3760                      }                      }
3761                      ///////////////                      ///////////////
3762                      // process y //                      // process y //
3763                      ///////////////                      ///////////////
3764                      if (add_EM_F) {                      if (add_EM_F) {
3765                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3766                          const double y_0 = y_p[0];                          const double tmp0_1 = w3*y_p[0];
                         const double tmp0_1 = w3*y_0;  
3767                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
3768                          EM_F[3]+=tmp0_1;                          EM_F[3]+=tmp0_1;
3769                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 BOTTOM */  
3770                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_N0*(m_N1-2)+k0;
3771                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3772                  }                  }
# Line 3810  void Rectangle::assemblePDEBoundarySyste Line 3788  void Rectangle::assemblePDEBoundarySyste
3788          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
3789          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
3790      }      }
     /* GENERATOR SNIP_PDEBC_SYSTEM_PRE TOP */  
3791      const double w0 = 0.31100423396407310779*h1;      const double w0 = 0.31100423396407310779*h1;
3792      const double w1 = 0.022329099369260225539*h1;      const double w1 = 0.022329099369260225539*h1;
3793      const double w10 = 0.022329099369260225539*h0;      const double w10 = 0.022329099369260225539*h0;
# Line 3827  void Rectangle::assemblePDEBoundarySyste Line 3804  void Rectangle::assemblePDEBoundarySyste
3804      const double w7 = 0.5*h1;      const double w7 = 0.5*h1;
3805      const double w8 = 0.083333333333333333333*h0;      const double w8 = 0.083333333333333333333*h0;
3806      const double w9 = 0.31100423396407310779*h0;      const double w9 = 0.31100423396407310779*h0;
     /* GENERATOR SNIP_PDEBC_SYSTEM_PRE BOTTOM */  
3807      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3808      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3809      rhs.requireWrite();      rhs.requireWrite();
# Line 3840  void Rectangle::assemblePDEBoundarySyste Line 3816  void Rectangle::assemblePDEBoundarySyste
3816                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
3817                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
3818                      const index_t e = k1;                      const index_t e = k1;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_0 TOP */  
3819                      ///////////////                      ///////////////
3820                      // process d //                      // process d //
3821                      ///////////////                      ///////////////
# Line 3852  void Rectangle::assemblePDEBoundarySyste Line 3827  void Rectangle::assemblePDEBoundarySyste
3827                                      const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
3828                                      const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
3829                                      const double tmp0_0 = d_0 + d_1;                                      const double tmp0_0 = d_0 + d_1;
                                     const double tmp1_1 = d_1*w1;  
                                     const double tmp4_1 = d_0*w1;  
3830                                      const double tmp0_1 = d_0*w0;                                      const double tmp0_1 = d_0*w0;
3831                                      const double tmp3_1 = d_1*w0;                                      const double tmp1_1 = d_1*w1;
3832                                      const double tmp2_1 = tmp0_0*w2;                                      const double tmp2_1 = tmp0_0*w2;
3833                                        const double tmp3_1 = d_1*w0;
3834                                        const double tmp4_1 = d_0*w1;
3835                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1;  
3836                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1;
3837                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1;
3838                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
3839                                  }                                  }
3840                               }                               }
# Line 3867  void Rectangle::assemblePDEBoundarySyste Line 3842  void Rectangle::assemblePDEBoundarySyste
3842                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3843                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3844                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
                                     const double tmp1_1 = d_0*w4;  
3845                                      const double tmp0_1 = d_0*w3;                                      const double tmp0_1 = d_0*w3;
3846                                        const double tmp1_1 = d_0*w4;
3847                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp1_1;  
3848                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1;
3849                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp1_1;
3850                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
3851                                  }                                  }
3852                              }                              }
# Line 3902  void Rectangle::assemblePDEBoundarySyste Line 3877  void Rectangle::assemblePDEBoundarySyste
3877                              }                              }
3878                          }                          }
3879                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_0 BOTTOM */  
3880                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_N0*k1;
3881                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
3882                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 3917  void Rectangle::assemblePDEBoundarySyste Line 3891  void Rectangle::assemblePDEBoundarySyste
3891                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
3892                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
3893                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_1 TOP */  
3894                      ///////////////                      ///////////////
3895                      // process d //                      // process d //
3896                      ///////////////                      ///////////////
# Line 3934  void Rectangle::assemblePDEBoundarySyste Line 3907  void Rectangle::assemblePDEBoundarySyste
3907                                      const double tmp3_1 = d_0*w0;                                      const double tmp3_1 = d_0*w0;
3908                                      const double tmp0_1 = d_1*w0;                                      const double tmp0_1 = d_1*w0;
3909                                      const double tmp2_1 = tmp0_0*w2;                                      const double tmp2_1 = tmp0_0*w2;
3910                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
3911                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1;
3912                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1;
3913                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3914                                  }                                  }
3915                               }                               }
3916                          } else { /* constant data */                          } else { /* constant data */
# Line 3946  void Rectangle::assemblePDEBoundarySyste Line 3919  void Rectangle::assemblePDEBoundarySyste
3919                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
3920                                      const double tmp1_1 = d_0*w4;                                      const double tmp1_1 = d_0*w4;
3921                                      const double tmp0_1 = d_0*w3;                                      const double tmp0_1 = d_0*w3;
3922                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
3923                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp1_1;
3924                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1;
3925                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
3926                                  }                                  }
3927                              }                              }
3928                          }                          }
# Line 3979  void Rectangle::assemblePDEBoundarySyste Line 3952  void Rectangle::assemblePDEBoundarySyste
3952                              }                              }
3953                          }                          }
3954                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_1 BOTTOM */  
3955                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_N0*(k1+1)-2;
3956                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
3957                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 3994  void Rectangle::assemblePDEBoundarySyste Line 3966  void Rectangle::assemblePDEBoundarySyste
3966                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
3967                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
3968                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_2 TOP */  
3969                      ///////////////                      ///////////////
3970                      // process d //                      // process d //
3971                      ///////////////                      ///////////////
# Line 4011  void Rectangle::assemblePDEBoundarySyste Line 3982  void Rectangle::assemblePDEBoundarySyste
3982                                      const double tmp0_1 = tmp0_0*w8;                                      const double tmp0_1 = tmp0_0*w8;
3983                                      const double tmp1_1 = d_1*w10;                                      const double tmp1_1 = d_1*w10;
3984                                      const double tmp3_1 = d_0*w10;                                      const double tmp3_1 = d_0*w10;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
3985                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3986                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
3987                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
3988                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
3989                                  }                                  }
3990                               }                               }
# Line 4023  void Rectangle::assemblePDEBoundarySyste Line 3994  void Rectangle::assemblePDEBoundarySyste
3994                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
3995                                      const double tmp0_1 = d_0*w11;                                      const double tmp0_1 = d_0*w11;
3996                                      const double tmp1_1 = d_0*w12;                                      const double tmp1_1 = d_0*w12;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
3997                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1;
3998                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
3999                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
4000                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1;
4001                                  }                                  }
4002                              }                              }
# Line 4056  void Rectangle::assemblePDEBoundarySyste Line 4027  void Rectangle::assemblePDEBoundarySyste
4027                              }                              }
4028                          }                          }
4029                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_2 BOTTOM */  
4030                      const index_t firstNode=k0;                      const index_t firstNode=k0;
4031                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4032                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 4071  void Rectangle::assemblePDEBoundarySyste Line 4041  void Rectangle::assemblePDEBoundarySyste
4041                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4042                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4043                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_3 TOP */  
4044                      ///////////////                      ///////////////
4045                      // process d //                      // process d //
4046                      ///////////////                      ///////////////
# Line 4088  void Rectangle::assemblePDEBoundarySyste Line 4057  void Rectangle::assemblePDEBoundarySyste
4057                                      const double tmp0_1 = tmp0_0*w8;                                      const double tmp0_1 = tmp0_0*w8;
4058                                      const double tmp3_1 = d_1*w10;                                      const double tmp3_1 = d_1*w10;
4059                                      const double tmp1_1 = d_0*w10;                                      const double tmp1_1 = d_0*w10;
4060                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
4061                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
4062                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
4063                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
4064                                  }                                  }
4065                               }                               }
4066                          } else { /* constant data */                          } else { /* constant data */
# Line 4100  void Rectangle::assemblePDEBoundarySyste Line 4069  void Rectangle::assemblePDEBoundarySyste
4069                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
4070                                      const double tmp0_1 = d_0*w11;                                      const double tmp0_1 = d_0*w11;
4071                                      const double tmp1_1 = d_0*w12;                                      const double tmp1_1 = d_0*w12;
4072                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp1_1;
4073                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1;  
4074                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
4075                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1;
4076                                  }                                  }
4077                              }                              }
4078                          }                          }
# Line 4133  void Rectangle::assemblePDEBoundarySyste Line 4102  void Rectangle::assemblePDEBoundarySyste
4102                              }                              }
4103                          }                          }
4104                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_3 BOTTOM */  
4105                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_N0*(m_N1-2)+k0;
4106                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4107                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 4156  void Rectangle::assemblePDEBoundarySyste Line 4124  void Rectangle::assemblePDEBoundarySyste
4124          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
4125          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
4126      }      }
     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_PRE TOP */  
4127      const double w0 = 0.25*h1;      const double w0 = 0.25*h1;
4128      const double w1 = 0.5*h1;      const double w1 = 0.5*h1;
4129      const double w2 = 0.25*h0;      const double w2 = 0.25*h0;
4130      const double w3 = 0.5*h0;      const double w3 = 0.5*h0;
     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_PRE BOTTOM */  
4131      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
4132      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
4133      rhs.requireWrite();      rhs.requireWrite();
# Line 4174  void Rectangle::assemblePDEBoundarySyste Line 4140  void Rectangle::assemblePDEBoundarySyste
4140                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4141                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4142                      const index_t e = k1;                      const index_t e = k1;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_0 TOP */  
4143                      ///////////////                      ///////////////
4144                      // process d //                      // process d //
4145                      ///////////////                      ///////////////
# Line 4185  void Rectangle::assemblePDEBoundarySyste Line 4150  void Rectangle::assemblePDEBoundarySyste
4150                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double d_0 = d_p[INDEX2(k, m, numEq)];
4151                                  const double tmp0_1 = d_0*w0;                                  const double tmp0_1 = d_0*w0;
4152                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;  
4153                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;
4154                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;
4155                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
4156                              }                              }
4157                          }                          }
# Line 4197  void Rectangle::assemblePDEBoundarySyste Line 4162  void Rectangle::assemblePDEBoundarySyste
4162                      if (add_EM_F) {                      if (add_EM_F) {
4163                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4164                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4165                              const double y_0 = y_p[k];                              const double tmp0_1 = w1*y_p[k];
                             const double tmp0_1 = w1*y_0;  
4166                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
4167                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
4168                          }                          }
4169                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_0 BOTTOM */  
4170                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_N0*k1;
4171                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4172                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
4173                  }                  }
                 // ADD EM_F  and EM_S  
4174              } // end colouring              } // end colouring
4175          }          }
4176    
# Line 4219  void Rectangle::assemblePDEBoundarySyste Line 4181  void Rectangle::assemblePDEBoundarySyste
4181                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4182                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4183                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_1 TOP */  
4184                      ///////////////                      ///////////////
4185                      // process d //                      // process d //
4186                      ///////////////                      ///////////////
# Line 4229  void Rectangle::assemblePDEBoundarySyste Line 4190  void Rectangle::assemblePDEBoundarySyste
4190                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
4191                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double d_0 = d_p[INDEX2(k, m, numEq)];
4192                                  const double tmp0_1 = d_0*w0;                                  const double tmp0_1 = d_0*w0;
4193                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
4194                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;
4195                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;
4196                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
4197                              }                              }
4198                          }                          }
4199                      }                      }
# Line 4242  void Rectangle::assemblePDEBoundarySyste Line 4203  void Rectangle::assemblePDEBoundarySyste
4203                      if (add_EM_F) {                      if (add_EM_F) {
4204                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4205                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4206                              const double y_0 = y_p[k];                              const double tmp0_1 = w1*y_p[k];
                             const double tmp0_1 = w1*y_0;  
4207                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
4208                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
4209                          }                          }
4210                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_1 BOTTOM */  
4211                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_N0*(k1+1)-2;
4212                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4213                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 4263  void Rectangle::assemblePDEBoundarySyste Line 4222  void Rectangle::assemblePDEBoundarySyste
4222                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4223                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4224                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_2 TOP */  
4225                      ///////////////                      ///////////////
4226                      // process d //                      // process d //
4227                      ///////////////                      ///////////////
# Line 4273  void Rectangle::assemblePDEBoundarySyste Line 4231  void Rectangle::assemblePDEBoundarySyste
4231                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
4232                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double d_0 = d_p[INDEX2(k, m, numEq)];
4233                                  const double tmp0_1 = d_0*w2;                                  const double tmp0_1 = d_0*w2;
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
4234                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
4235                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
4236                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
4237                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
4238                              }                              }
4239                          }                          }
# Line 4286  void Rectangle::assemblePDEBoundarySyste Line 4244  void Rectangle::assemblePDEBoundarySyste
4244                      if (add_EM_F) {                      if (add_EM_F) {
4245                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4246                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4247                              const double y_0 = y_p[k];                              const double tmp0_1 = w3*y_p[k];
                             const double tmp0_1 = w3*y_0;  
4248                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
4249                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
4250                          }                          }
4251                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_2 BOTTOM */  
4252                      const index_t firstNode=k0;                      const index_t firstNode=k0;
4253                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4254                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 4307  void Rectangle::assemblePDEBoundarySyste Line 4263  void Rectangle::assemblePDEBoundarySyste
4263                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4264                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4265                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_3 TOP */  
4266                      ///////////////                      ///////////////
4267                      // process d //                      // process d //
4268                      ///////////////                      ///////////////
# Line 4317  void Rectangle::assemblePDEBoundarySyste Line 4272  void Rectangle::assemblePDEBoundarySyste
4272                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
4273                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double d_0 = d_p[INDEX2(k, m, numEq)];
4274                                  const double tmp0_1 = d_0*w2;                                  const double tmp0_1 = d_0*w2;
4275                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
4276                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;  
4277                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
4278                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
4279                              }                              }
4280                          }                          }
4281                      }                      }
# Line 4330  void Rectangle::assemblePDEBoundarySyste Line 4285  void Rectangle::assemblePDEBoundarySyste
4285                      if (add_EM_F) {                      if (add_EM_F) {
4286                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4287                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4288                              const double y_0 = y_p[k];                              const double tmp0_1 = w3*y_p[k];
                             const double tmp0_1 = w3*y_0;  
4289                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
4290                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
4291                          }                          }
4292                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_3 BOTTOM */  
4293                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_N0*(m_N1-2)+k0;
4294                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4295                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);

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

  ViewVC Help
Powered by ViewVC 1.1.26