/[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 3769 by caltinay, Tue Jan 17 04:09:55 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 488  void Rectangle::setToSize(escript::Data& Line 492  void Rectangle::setToSize(escript::Data&
492  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
493                                                  bool reducedColOrder) const                                                  bool reducedColOrder) const
494  {  {
495        /* FIXME: reduced
496      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
497          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
498        */
499      return m_pattern;      return m_pattern;
500  }  }
501    
# Line 556  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 1272  void Rectangle::createPattern() Line 1277  void Rectangle::createPattern()
1277      Paso_Pattern_free(rowPattern);      Paso_Pattern_free(rowPattern);
1278  }  }
1279    
1280    //private
1281    void Rectangle::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
1282             const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
1283             bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
1284    {
1285        IndexVector rowIndex;
1286        rowIndex.push_back(m_dofMap[firstNode]);
1287        rowIndex.push_back(m_dofMap[firstNode+1]);
1288        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
1289        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
1290        if (addF) {
1291            double *F_p=F.getSampleDataRW(0);
1292            for (index_t i=0; i<rowIndex.size(); i++) {
1293                if (rowIndex[i]<getNumDOF()) {
1294                    for (index_t eq=0; eq<nEq; eq++) {
1295                        F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
1296                    }
1297                }
1298            }
1299        }
1300        if (addS) {
1301            addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
1302        }
1303    }
1304    
1305  //protected  //protected
1306  void Rectangle::interpolateNodesOnElements(escript::Data& out,  void Rectangle::interpolateNodesOnElements(escript::Data& out,
1307                                          escript::Data& in, bool reduced) const                                          escript::Data& in, bool reduced) const
# Line 2133  void Rectangle::assemblePDESingle(Paso_S Line 2163  void Rectangle::assemblePDESingle(Paso_S
2163    
2164                      // 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)
2165                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
2166                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
                     rowIndex.push_back(m_dofMap[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]);  
                     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);  
                     }  
   
2167                  } // end k0 loop                  } // end k0 loop
2168              } // end k1 loop              } // end k1 loop
2169          } // end of colouring          } // end of colouring
# Line 2339  void Rectangle::assemblePDESingleReduced Line 2353  void Rectangle::assemblePDESingleReduced
2353    
2354                      // 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)
2355                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
2356                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
                     rowIndex.push_back(m_dofMap[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]);  
                     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);  
                     }  
   
2357                  } // end k0 loop                  } // end k0 loop
2358              } // end k1 loop              } // end k1 loop
2359          } // end of colouring          } // end of colouring
# Line 2731  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 2762  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 3116  void Rectangle::assemblePDESystem(Paso_S Line 3114  void Rectangle::assemblePDESystem(Paso_S
3114    
3115                      // 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)
3116                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
3117                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
3118                      rowIndex.push_back(m_dofMap[firstNode]);                              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]);  
                     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);  
                     }  
   
3119                  } // end k0 loop                  } // end k0 loop
3120              } // end k1 loop              } // end k1 loop
3121          } // end of colouring          } // end of colouring
# Line 3358  void Rectangle::assemblePDESystemReduced Line 3339  void Rectangle::assemblePDESystemReduced
3339    
3340                      // 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)
3341                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
3342                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
3343                      rowIndex.push_back(m_dofMap[firstNode]);                              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]);  
                     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);  
                     }  
   
3344                  } // end k0 loop                  } // end k0 loop
3345              } // end k1 loop              } // end k1 loop
3346          } // end of colouring          } // end of colouring
# Line 3387  void Rectangle::assemblePDESystemReduced Line 3351  void Rectangle::assemblePDESystemReduced
3351  void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
3352        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3353  {  {
3354        const double h0 = m_l0/m_gNE0;
3355        const double h1 = m_l1/m_gNE1;
3356        const double w0 = 0.31100423396407310779*h1;
3357        const double w1 = 0.022329099369260225539*h1;
3358        const double w10 = 0.022329099369260225539*h0;
3359        const double w11 = 0.16666666666666666667*h0;
3360        const double w12 = 0.33333333333333333333*h0;
3361        const double w13 = 0.39433756729740644113*h0;
3362        const double w14 = 0.10566243270259355887*h0;
3363        const double w15 = 0.5*h0;
3364        const double w2 = 0.083333333333333333333*h1;
3365        const double w3 = 0.33333333333333333333*h1;
3366        const double w4 = 0.16666666666666666667*h1;
3367        const double w5 = 0.39433756729740644113*h1;
3368        const double w6 = 0.10566243270259355887*h1;
3369        const double w7 = 0.5*h1;
3370        const double w8 = 0.083333333333333333333*h0;
3371        const double w9 = 0.31100423396407310779*h0;
3372        const bool add_EM_S=!d.isEmpty();
3373        const bool add_EM_F=!y.isEmpty();
3374        rhs.requireWrite();
3375    #pragma omp parallel
3376        {
3377            if (m_faceOffset[0] > -1) {
3378                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3379    #pragma omp for nowait
3380                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3381                        vector<double> EM_S(4*4, 0);
3382                        vector<double> EM_F(4, 0);
3383                        const index_t e = k1;
3384                        ///////////////
3385                        // process d //
3386                        ///////////////
3387                        if (add_EM_S) {
3388                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3389                            if (d.actsExpanded()) {
3390                                const double d_0 = d_p[0];
3391                                const double d_1 = d_p[1];
3392                                const double tmp0_0 = d_0 + d_1;
3393                                const double tmp1_1 = d_1*w1;
3394                                const double tmp4_1 = d_0*w1;
3395                                const double tmp0_1 = d_0*w0;
3396                                const double tmp3_1 = d_1*w0;
3397                                const double tmp2_1 = tmp0_0*w2;
3398                                EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
3399                                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;
3402                            } else { /* constant data */
3403                                const double d_0 = d_p[0];
3404                                const double tmp1_1 = d_0*w4;
3405                                const double tmp0_1 = d_0*w3;
3406                                EM_S[INDEX2(0,0,4)]+=tmp0_1;
3407                                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;
3410                            }
3411                        }
3412                        ///////////////
3413                        // process y //
3414                        ///////////////
3415                        if (add_EM_F) {
3416                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3417                            if (y.actsExpanded()) {
3418                                const double y_0 = y_p[0];
3419                                const double y_1 = y_p[1];
3420                                const double tmp3_1 = w5*y_1;
3421                                const double tmp2_1 = w6*y_0;
3422                                const double tmp0_1 = w6*y_1;
3423                                const double tmp1_1 = w5*y_0;
3424                                EM_F[0]+=tmp0_1 + tmp1_1;
3425                                EM_F[2]+=tmp2_1 + tmp3_1;
3426                            } else { /* constant data */
3427                                const double y_0 = y_p[0];
3428                                const double tmp0_1 = w7*y_0;
3429                                EM_F[0]+=tmp0_1;
3430                                EM_F[2]+=tmp0_1;
3431                            }
3432                        }
3433                        const index_t firstNode=m_N0*k1;
3434                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3435                    }
3436                } // end colouring
3437            }
3438    
3439            if (m_faceOffset[1] > -1) {
3440                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3441    #pragma omp for nowait
3442                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3443                        vector<double> EM_S(4*4, 0);
3444                        vector<double> EM_F(4, 0);
3445                        const index_t e = m_faceOffset[1]+k1;
3446                        ///////////////
3447                        // process d //
3448                        ///////////////
3449                        if (add_EM_S) {
3450                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3451                            if (d.actsExpanded()) {
3452                                const double d_0 = d_p[0];
3453                                const double d_1 = d_p[1];
3454                                const double tmp0_0 = d_0 + d_1;
3455                                const double tmp4_1 = d_1*w1;
3456                                const double tmp1_1 = d_0*w1;
3457                                const double tmp3_1 = d_0*w0;
3458                                const double tmp0_1 = d_1*w0;
3459                                const double tmp2_1 = tmp0_0*w2;
3460                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;
3461                                EM_S[INDEX2(3,1,4)]+=tmp2_1;
3462                                EM_S[INDEX2(1,3,4)]+=tmp2_1;
3463                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp1_1;
3464                            } else { /* constant data */
3465                                const double d_0 = d_p[0];
3466                                const double tmp1_1 = d_0*w4;
3467                                const double tmp0_1 = d_0*w3;
3468                                EM_S[INDEX2(1,1,4)]+=tmp0_1;
3469                                EM_S[INDEX2(3,1,4)]+=tmp1_1;
3470                                EM_S[INDEX2(1,3,4)]+=tmp1_1;
3471                                EM_S[INDEX2(3,3,4)]+=tmp0_1;
3472                            }
3473                        }
3474                        ///////////////
3475                        // process y //
3476                        ///////////////
3477                        if (add_EM_F) {
3478                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3479                            if (y.actsExpanded()) {
3480                                const double y_0 = y_p[0];
3481                                const double y_1 = y_p[1];
3482                                const double tmp3_1 = w5*y_1;
3483                                const double tmp2_1 = w6*y_0;
3484                                const double tmp0_1 = w6*y_1;
3485                                const double tmp1_1 = w5*y_0;
3486                                EM_F[1]+=tmp0_1 + tmp1_1;
3487                                EM_F[3]+=tmp2_1 + tmp3_1;
3488                            } else { /* constant data */
3489                                const double y_0 = y_p[0];
3490                                const double tmp0_1 = w7*y_0;
3491                                EM_F[1]+=tmp0_1;
3492                                EM_F[3]+=tmp0_1;
3493                            }
3494                        }
3495                        const index_t firstNode=m_N0*(k1+1)-2;
3496                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3497                    }
3498                } // end colouring
3499            }
3500    
3501            if (m_faceOffset[2] > -1) {
3502                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3503    #pragma omp for nowait
3504                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3505                        vector<double> EM_S(4*4, 0);
3506                        vector<double> EM_F(4, 0);
3507                        const index_t e = m_faceOffset[2]+k0;
3508                        ///////////////
3509                        // process d //
3510                        ///////////////
3511                        if (add_EM_S) {
3512                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3513                            if (d.actsExpanded()) {
3514                                const double d_0 = d_p[0];
3515                                const double d_1 = d_p[1];
3516                                const double tmp0_0 = d_0 + d_1;
3517                                const double tmp4_1 = d_1*w9;
3518                                const double tmp2_1 = d_0*w9;
3519                                const double tmp0_1 = tmp0_0*w8;
3520                                const double tmp1_1 = d_1*w10;
3521                                const double tmp3_1 = d_0*w10;
3522                                EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;
3523                                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;
3526                            } else { /* constant data */
3527                                const double d_0 = d_p[0];
3528                                const double tmp0_1 = d_0*w11;
3529                                const double tmp1_1 = d_0*w12;
3530                                EM_S[INDEX2(0,0,4)]+=tmp1_1;
3531                                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;
3534                            }
3535                        }
3536                        ///////////////
3537                        // process y //
3538                        ///////////////
3539                        if (add_EM_F) {
3540                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3541                            if (y.actsExpanded()) {
3542                                const double y_0 = y_p[0];
3543                                const double y_1 = y_p[1];
3544                                const double tmp2_1 = w13*y_1;
3545                                const double tmp1_1 = w14*y_1;
3546                                const double tmp3_1 = w14*y_0;
3547                                const double tmp0_1 = w13*y_0;
3548                                EM_F[0]+=tmp0_1 + tmp1_1;
3549                                EM_F[1]+=tmp2_1 + tmp3_1;
3550                            } else { /* constant data */
3551                                const double y_0 = y_p[0];
3552                                const double tmp0_1 = w15*y_0;
3553                                EM_F[0]+=tmp0_1;
3554                                EM_F[1]+=tmp0_1;
3555                            }
3556                        }
3557                        const index_t firstNode=k0;
3558                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3559                    }
3560                } // end colouring
3561            }
3562    
3563            if (m_faceOffset[3] > -1) {
3564                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3565    #pragma omp for nowait
3566                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3567                        const index_t e = m_faceOffset[3]+k0;
3568                        vector<double> EM_S(4*4, 0);
3569                        vector<double> EM_F(4, 0);
3570                        ///////////////
3571                        // process d //
3572                        ///////////////
3573                        if (add_EM_S) {
3574                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3575                            if (d.actsExpanded()) {
3576                                const double d_0 = d_p[0];
3577                                const double d_1 = d_p[1];
3578                                const double tmp0_0 = d_0 + d_1;
3579                                const double tmp2_1 = d_1*w9;
3580                                const double tmp4_1 = d_0*w9;
3581                                const double tmp0_1 = tmp0_0*w8;
3582                                const double tmp3_1 = d_1*w10;
3583                                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;
3586                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
3587                                EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
3588                            } else { /* constant data */
3589                                const double d_0 = d_p[0];
3590                                const double tmp0_1 = d_0*w11;
3591                                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;
3594                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
3595                                EM_S[INDEX2(3,3,4)]+=tmp1_1;
3596                            }
3597                        }
3598                        ///////////////
3599                        // process y //
3600                        ///////////////
3601                        if (add_EM_F) {
3602                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3603                            if (y.actsExpanded()) {
3604                                const double y_0 = y_p[0];
3605                                const double y_1 = y_p[1];
3606                                const double tmp2_1 = w13*y_1;
3607                                const double tmp1_1 = w14*y_1;
3608                                const double tmp3_1 = w14*y_0;
3609                                const double tmp0_1 = w13*y_0;
3610                                EM_F[2]+=tmp0_1 + tmp1_1;
3611                                EM_F[3]+=tmp2_1 + tmp3_1;
3612                            } else { /* constant data */
3613                                const double y_0 = y_p[0];
3614                                const double tmp0_1 = w15*y_0;
3615                                EM_F[2]+=tmp0_1;
3616                                EM_F[3]+=tmp0_1;
3617                            }
3618                        }
3619                        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);
3621                    }
3622                } // end colouring
3623            }
3624        } // end of parallel section
3625  }  }
3626    
3627  //protected  //protected
3628  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
3629        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3630  {  {
3631        const double h0 = m_l0/m_gNE0;
3632        const double h1 = m_l1/m_gNE1;
3633        const double w0 = 0.25*h1;
3634        const double w1 = 0.5*h1;
3635        const double w2 = 0.25*h0;
3636        const double w3 = 0.5*h0;
3637        const bool add_EM_S=!d.isEmpty();
3638        const bool add_EM_F=!y.isEmpty();
3639        rhs.requireWrite();
3640    #pragma omp parallel
3641        {
3642            if (m_faceOffset[0] > -1) {
3643                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3644    #pragma omp for nowait
3645                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3646                        vector<double> EM_S(4*4, 0);
3647                        vector<double> EM_F(4, 0);
3648                        const index_t e = k1;
3649                        ///////////////
3650                        // process d //
3651                        ///////////////
3652                        if (add_EM_S) {
3653                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3654                            const double d_0 = d_p[0];
3655                            const double tmp0_1 = d_0*w0;
3656                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
3657                            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;
3660                        }
3661                        ///////////////
3662                        // process y //
3663                        ///////////////
3664                        if (add_EM_F) {
3665                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3666                            const double tmp0_1 = w1*y_p[0];
3667                            EM_F[0]+=tmp0_1;
3668                            EM_F[2]+=tmp0_1;
3669                        }
3670                        const index_t firstNode=m_N0*k1;
3671                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3672                    }
3673                } // end colouring
3674            }
3675    
3676            if (m_faceOffset[1] > -1) {
3677                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3678    #pragma omp for nowait
3679                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3680                        vector<double> EM_S(4*4, 0);
3681                        vector<double> EM_F(4, 0);
3682                        const index_t e = m_faceOffset[1]+k1;
3683                        ///////////////
3684                        // process d //
3685                        ///////////////
3686                        if (add_EM_S) {
3687                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3688                            const double d_0 = d_p[0];
3689                            const double tmp0_1 = d_0*w0;
3690                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
3691                            EM_S[INDEX2(3,1,4)]+=tmp0_1;
3692                            EM_S[INDEX2(1,3,4)]+=tmp0_1;
3693                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
3694                        }
3695                        ///////////////
3696                        // process y //
3697                        ///////////////
3698                        if (add_EM_F) {
3699                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3700                            const double tmp0_1 = w1*y_p[0];
3701                            EM_F[1]+=tmp0_1;
3702                            EM_F[3]+=tmp0_1;
3703                        }
3704                        const index_t firstNode=m_N0*(k1+1)-2;
3705                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3706                    }
3707                } // end colouring
3708            }
3709    
3710            if (m_faceOffset[2] > -1) {
3711                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3712    #pragma omp for nowait
3713                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3714                        vector<double> EM_S(4*4, 0);
3715                        vector<double> EM_F(4, 0);
3716                        const index_t e = m_faceOffset[2]+k0;
3717                        ///////////////
3718                        // process d //
3719                        ///////////////
3720                        if (add_EM_S) {
3721                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3722                            const double tmp0_1 = d_p[0]*w2;
3723                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
3724                            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;
3727                        }
3728                        ///////////////
3729                        // process y //
3730                        ///////////////
3731                        if (add_EM_F) {
3732                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3733                            const double tmp0_1 = w3*y_p[0];
3734                            EM_F[0]+=tmp0_1;
3735                            EM_F[1]+=tmp0_1;
3736                        }
3737                        const index_t firstNode=k0;
3738                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3739                    }
3740                } // end colouring
3741            }
3742    
3743            if (m_faceOffset[3] > -1) {
3744                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3745    #pragma omp for nowait
3746                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3747                        vector<double> EM_S(4*4, 0);
3748                        vector<double> EM_F(4, 0);
3749                        const index_t e = m_faceOffset[3]+k0;
3750                        ///////////////
3751                        // process d //
3752                        ///////////////
3753                        if (add_EM_S) {
3754                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3755                            const double tmp0_1 = d_p[0]*w2;
3756                            EM_S[INDEX2(2,2,4)]+=tmp0_1;
3757                            EM_S[INDEX2(3,2,4)]+=tmp0_1;
3758                            EM_S[INDEX2(2,3,4)]+=tmp0_1;
3759                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
3760                        }
3761                        ///////////////
3762                        // process y //
3763                        ///////////////
3764                        if (add_EM_F) {
3765                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3766                            const double tmp0_1 = w3*y_p[0];
3767                            EM_F[2]+=tmp0_1;
3768                            EM_F[3]+=tmp0_1;
3769                        }
3770                        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);
3772                    }
3773                } // end colouring
3774            }
3775        } // end of parallel section
3776  }  }
3777    
3778  //protected  //protected
3779  void Rectangle::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
3780        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3781  {  {
3782        const double h0 = m_l0/m_gNE0;
3783        const double h1 = m_l1/m_gNE1;
3784        dim_t numEq, numComp;
3785        if (!mat) {
3786            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
3787        } else {
3788            numEq=mat->logical_row_block_size;
3789            numComp=mat->logical_col_block_size;
3790        }
3791        const double w0 = 0.31100423396407310779*h1;
3792        const double w1 = 0.022329099369260225539*h1;
3793        const double w10 = 0.022329099369260225539*h0;
3794        const double w11 = 0.16666666666666666667*h0;
3795        const double w12 = 0.33333333333333333333*h0;
3796        const double w13 = 0.39433756729740644113*h0;
3797        const double w14 = 0.10566243270259355887*h0;
3798        const double w15 = 0.5*h0;
3799        const double w2 = 0.083333333333333333333*h1;
3800        const double w3 = 0.33333333333333333333*h1;
3801        const double w4 = 0.16666666666666666667*h1;
3802        const double w5 = 0.39433756729740644113*h1;
3803        const double w6 = 0.10566243270259355887*h1;
3804        const double w7 = 0.5*h1;
3805        const double w8 = 0.083333333333333333333*h0;
3806        const double w9 = 0.31100423396407310779*h0;
3807        const bool add_EM_S=!d.isEmpty();
3808        const bool add_EM_F=!y.isEmpty();
3809        rhs.requireWrite();
3810    #pragma omp parallel
3811        {
3812            if (m_faceOffset[0] > -1) {
3813                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3814    #pragma omp for nowait
3815                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3816                        vector<double> EM_S(4*4*numEq*numComp, 0);
3817                        vector<double> EM_F(4*numEq, 0);
3818                        const index_t e = k1;
3819                        ///////////////
3820                        // process d //
3821                        ///////////////
3822                        if (add_EM_S) {
3823                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3824                            if (d.actsExpanded()) {
3825                                for (index_t k=0; k<numEq; k++) {
3826                                    for (index_t m=0; m<numComp; m++) {
3827                                        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)];
3829                                        const double tmp0_0 = d_0 + d_1;
3830                                        const double tmp0_1 = d_0*w0;
3831                                        const double tmp1_1 = d_1*w1;
3832                                        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;
3836                                        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;
3839                                    }
3840                                 }
3841                            } else { /* constant data */
3842                                for (index_t k=0; k<numEq; k++) {
3843                                    for (index_t m=0; m<numComp; m++) {
3844                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
3845                                        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;
3848                                        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;
3851                                    }
3852                                }
3853                            }
3854                        }
3855                        ///////////////
3856                        // process y //
3857                        ///////////////
3858                        if (add_EM_F) {
3859                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3860                            if (y.actsExpanded()) {
3861                                for (index_t k=0; k<numEq; k++) {
3862                                    const double y_0 = y_p[INDEX2(k, 0, numEq)];
3863                                    const double y_1 = y_p[INDEX2(k, 1, numEq)];
3864                                    const double tmp3_1 = w5*y_1;
3865                                    const double tmp2_1 = w6*y_0;
3866                                    const double tmp0_1 = w6*y_1;
3867                                    const double tmp1_1 = w5*y_0;
3868                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3869                                    EM_F[INDEX2(k,2,numEq)]+=tmp2_1 + tmp3_1;
3870                                }
3871                            } else { /* constant data */
3872                                for (index_t k=0; k<numEq; k++) {
3873                                    const double y_0 = y_p[k];
3874                                    const double tmp0_1 = w7*y_0;
3875                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3876                                    EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
3877                                }
3878                            }
3879                        }
3880                        const index_t firstNode=m_N0*k1;
3881                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
3882                                firstNode, numEq, numComp);
3883                    }
3884                } // end colouring
3885            }
3886    
3887            if (m_faceOffset[1] > -1) {
3888                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3889    #pragma omp for nowait
3890                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3891                        vector<double> EM_S(4*4*numEq*numComp, 0);
3892                        vector<double> EM_F(4*numEq, 0);
3893                        const index_t e = m_faceOffset[1]+k1;
3894                        ///////////////
3895                        // process d //
3896                        ///////////////
3897                        if (add_EM_S) {
3898                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3899                            if (d.actsExpanded()) {
3900                                for (index_t k=0; k<numEq; k++) {
3901                                    for (index_t m=0; m<numComp; m++) {
3902                                        const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
3903                                        const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
3904                                        const double tmp0_0 = d_0 + d_1;
3905                                        const double tmp4_1 = d_1*w1;
3906                                        const double tmp1_1 = d_0*w1;
3907                                        const double tmp3_1 = d_0*w0;
3908                                        const double tmp0_1 = d_1*w0;
3909                                        const double tmp2_1 = tmp0_0*w2;
3910                                        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;
3912                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1;
3913                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3914                                    }
3915                                 }
3916                            } else { /* constant data */
3917                                for (index_t k=0; k<numEq; k++) {
3918                                    for (index_t m=0; m<numComp; m++) {
3919                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
3920                                        const double tmp1_1 = d_0*w4;
3921                                        const double tmp0_1 = d_0*w3;
3922                                        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;
3924                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1;
3925                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
3926                                    }
3927                                }
3928                            }
3929                        }
3930                        ///////////////
3931                        // process y //
3932                        ///////////////
3933                        if (add_EM_F) {
3934                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3935                            if (y.actsExpanded()) {
3936                                for (index_t k=0; k<numEq; k++) {
3937                                    const double y_0 = y_p[INDEX2(k, 0, numEq)];
3938                                    const double y_1 = y_p[INDEX2(k, 1, numEq)];
3939                                    const double tmp3_1 = w5*y_1;
3940                                    const double tmp2_1 = w6*y_0;
3941                                    const double tmp0_1 = w6*y_1;
3942                                    const double tmp1_1 = w5*y_0;
3943                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp1_1;
3944                                    EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;
3945                                }
3946                            } else { /* constant data */
3947                                for (index_t k=0; k<numEq; k++) {
3948                                    const double y_0 = y_p[k];
3949                                    const double tmp0_1 = w7*y_0;
3950                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3951                                    EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
3952                                }
3953                            }
3954                        }
3955                        const index_t firstNode=m_N0*(k1+1)-2;
3956                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
3957                                firstNode, numEq, numComp);
3958                    }
3959                } // end colouring
3960            }
3961    
3962            if (m_faceOffset[2] > -1) {
3963                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3964    #pragma omp for nowait
3965                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3966                        vector<double> EM_S(4*4*numEq*numComp, 0);
3967                        vector<double> EM_F(4*numEq, 0);
3968                        const index_t e = m_faceOffset[2]+k0;
3969                        ///////////////
3970                        // process d //
3971                        ///////////////
3972                        if (add_EM_S) {
3973                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3974                            if (d.actsExpanded()) {
3975                                for (index_t k=0; k<numEq; k++) {
3976                                    for (index_t m=0; m<numComp; m++) {
3977                                        const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
3978                                        const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
3979                                        const double tmp0_0 = d_0 + d_1;
3980                                        const double tmp4_1 = d_1*w9;
3981                                        const double tmp2_1 = d_0*w9;
3982                                        const double tmp0_1 = tmp0_0*w8;
3983                                        const double tmp1_1 = d_1*w10;
3984                                        const double tmp3_1 = d_0*w10;
3985                                        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;
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;
3989                                    }
3990                                 }
3991                            } else { /* constant data */
3992                                for (index_t k=0; k<numEq; k++) {
3993                                    for (index_t m=0; m<numComp; m++) {
3994                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
3995                                        const double tmp0_1 = d_0*w11;
3996                                        const double tmp1_1 = d_0*w12;
3997                                        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;
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;
4001                                    }
4002                                }
4003                            }
4004                        }
4005                        ///////////////
4006                        // process y //
4007                        ///////////////
4008                        if (add_EM_F) {
4009                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4010                            if (y.actsExpanded()) {
4011                                for (index_t k=0; k<numEq; k++) {
4012                                    const double y_0 = y_p[INDEX2(k, 0, numEq)];
4013                                    const double y_1 = y_p[INDEX2(k, 1, numEq)];
4014                                    const double tmp2_1 = w13*y_1;
4015                                    const double tmp1_1 = w14*y_1;
4016                                    const double tmp3_1 = w14*y_0;
4017                                    const double tmp0_1 = w13*y_0;
4018                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
4019                                    EM_F[INDEX2(k,1,numEq)]+=tmp2_1 + tmp3_1;
4020                                }
4021                            } else { /* constant data */
4022                                for (index_t k=0; k<numEq; k++) {
4023                                    const double y_0 = y_p[k];
4024                                    const double tmp0_1 = w15*y_0;
4025                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
4026                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
4027                                }
4028                            }
4029                        }
4030                        const index_t firstNode=k0;
4031                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4032                                firstNode, numEq, numComp);
4033                    }
4034                } // end colouring
4035            }
4036    
4037            if (m_faceOffset[3] > -1) {
4038                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4039    #pragma omp for nowait
4040                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4041                        vector<double> EM_S(4*4*numEq*numComp, 0);
4042                        vector<double> EM_F(4*numEq, 0);
4043                        const index_t e = m_faceOffset[3]+k0;
4044                        ///////////////
4045                        // process d //
4046                        ///////////////
4047                        if (add_EM_S) {
4048                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
4049                            if (d.actsExpanded()) {
4050                                for (index_t k=0; k<numEq; k++) {
4051                                    for (index_t m=0; m<numComp; m++) {
4052                                        const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
4053                                        const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
4054                                        const double tmp0_0 = d_0 + d_1;
4055                                        const double tmp2_1 = d_1*w9;
4056                                        const double tmp4_1 = d_0*w9;
4057                                        const double tmp0_1 = tmp0_0*w8;
4058                                        const double tmp3_1 = d_1*w10;
4059                                        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;
4062                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
4063                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
4064                                    }
4065                                 }
4066                            } else { /* constant data */
4067                                for (index_t k=0; k<numEq; k++) {
4068                                    for (index_t m=0; m<numComp; m++) {
4069                                        const double d_0 = d_p[INDEX2(k, m, numEq)];
4070                                        const double tmp0_1 = d_0*w11;
4071                                        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;
4074                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
4075                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1;
4076                                    }
4077                                }
4078                            }
4079                        }
4080                        ///////////////
4081                        // process y //
4082                        ///////////////
4083                        if (add_EM_F) {
4084                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4085                            if (y.actsExpanded()) {
4086                                for (index_t k=0; k<numEq; k++) {
4087                                    const double y_0 = y_p[INDEX2(k, 0, numEq)];
4088                                    const double y_1 = y_p[INDEX2(k, 1, numEq)];
4089                                    const double tmp2_1 = w13*y_1;
4090                                    const double tmp1_1 = w14*y_1;
4091                                    const double tmp3_1 = w14*y_0;
4092                                    const double tmp0_1 = w13*y_0;
4093                                    EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp1_1;
4094                                    EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;
4095                                }
4096                            } else { /* constant data */
4097                                for (index_t k=0; k<numEq; k++) {
4098                                    const double y_0 = y_p[k];
4099                                    const double tmp0_1 = w15*y_0;
4100                                    EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
4101                                    EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
4102                                }
4103                            }
4104                        }
4105                        const index_t firstNode=m_N0*(m_N1-2)+k0;
4106                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4107                                firstNode, numEq, numComp);
4108                    }
4109                } // end colouring
4110            }
4111        } // end of parallel section
4112  }  }
4113    
4114  //protected  //protected
4115  void Rectangle::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
4116        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
4117  {  {
4118        const double h0 = m_l0/m_gNE0;
4119        const double h1 = m_l1/m_gNE1;
4120        dim_t numEq, numComp;
4121        if (!mat)
4122            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
4123        else {
4124            numEq=mat->logical_row_block_size;
4125            numComp=mat->logical_col_block_size;
4126        }
4127        const double w0 = 0.25*h1;
4128        const double w1 = 0.5*h1;
4129        const double w2 = 0.25*h0;
4130        const double w3 = 0.5*h0;
4131        const bool add_EM_S=!d.isEmpty();
4132        const bool add_EM_F=!y.isEmpty();
4133        rhs.requireWrite();
4134    #pragma omp parallel
4135        {
4136            if (m_faceOffset[0] > -1) {
4137                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
4138    #pragma omp for nowait
4139                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4140                        vector<double> EM_S(4*4*numEq*numComp, 0);
4141                        vector<double> EM_F(4*numEq, 0);
4142                        const index_t e = k1;
4143                        ///////////////
4144                        // process d //
4145                        ///////////////
4146                        if (add_EM_S) {
4147                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
4148                            for (index_t k=0; k<numEq; k++) {
4149                                for (index_t m=0; m<numComp; m++) {
4150                                    const double d_0 = d_p[INDEX2(k, m, numEq)];
4151                                    const double tmp0_1 = d_0*w0;
4152                                    EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
4153                                    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;
4156                                }
4157                            }
4158                        }
4159                        ///////////////
4160                        // process y //
4161                        ///////////////
4162                        if (add_EM_F) {
4163                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4164                            for (index_t k=0; k<numEq; k++) {
4165                                const double tmp0_1 = w1*y_p[k];
4166                                EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
4167                                EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
4168                            }
4169                        }
4170                        const index_t firstNode=m_N0*k1;
4171                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4172                                firstNode, numEq, numComp);
4173                    }
4174                } // end colouring
4175            }
4176    
4177            if (m_faceOffset[1] > -1) {
4178                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4179    #pragma omp for nowait
4180                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4181                        vector<double> EM_S(4*4*numEq*numComp, 0);
4182                        vector<double> EM_F(4*numEq, 0);
4183                        const index_t e = m_faceOffset[1]+k1;
4184                        ///////////////
4185                        // process d //
4186                        ///////////////
4187                        if (add_EM_S) {
4188                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
4189                            for (index_t k=0; k<numEq; k++) {
4190                                for (index_t m=0; m<numComp; m++) {
4191                                    const double d_0 = d_p[INDEX2(k, m, numEq)];
4192                                    const double tmp0_1 = d_0*w0;
4193                                    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;
4195                                    EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;
4196                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
4197                                }
4198                            }
4199                        }
4200                        ///////////////
4201                        // process y //
4202                        ///////////////
4203                        if (add_EM_F) {
4204                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4205                            for (index_t k=0; k<numEq; k++) {
4206                                const double tmp0_1 = w1*y_p[k];
4207                                EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
4208                                EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
4209                            }
4210                        }
4211                        const index_t firstNode=m_N0*(k1+1)-2;
4212                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4213                                firstNode, numEq, numComp);
4214                    }
4215                } // end colouring
4216            }
4217    
4218            if (m_faceOffset[2] > -1) {
4219                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4220    #pragma omp for nowait
4221                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4222                        vector<double> EM_S(4*4*numEq*numComp, 0);
4223                        vector<double> EM_F(4*numEq, 0);
4224                        const index_t e = m_faceOffset[2]+k0;
4225                        ///////////////
4226                        // process d //
4227                        ///////////////
4228                        if (add_EM_S) {
4229                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
4230                            for (index_t k=0; k<numEq; k++) {
4231                                for (index_t m=0; m<numComp; m++) {
4232                                    const double d_0 = d_p[INDEX2(k, m, numEq)];
4233                                    const double tmp0_1 = d_0*w2;
4234                                    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;
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;
4238                                }
4239                            }
4240                        }
4241                        ///////////////
4242                        // process y //
4243                        ///////////////
4244                        if (add_EM_F) {
4245                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4246                            for (index_t k=0; k<numEq; k++) {
4247                                const double tmp0_1 = w3*y_p[k];
4248                                EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
4249                                EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
4250                            }
4251                        }
4252                        const index_t firstNode=k0;
4253                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4254                                firstNode, numEq, numComp);
4255                    }
4256                } // end colouring
4257            }
4258    
4259            if (m_faceOffset[3] > -1) {
4260                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4261    #pragma omp for nowait
4262                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4263                        vector<double> EM_S(4*4*numEq*numComp, 0);
4264                        vector<double> EM_F(4*numEq, 0);
4265                        const index_t e = m_faceOffset[3]+k0;
4266                        ///////////////
4267                        // process d //
4268                        ///////////////
4269                        if (add_EM_S) {
4270                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
4271                            for (index_t k=0; k<numEq; k++) {
4272                                for (index_t m=0; m<numComp; m++) {
4273                                    const double d_0 = d_p[INDEX2(k, m, numEq)];
4274                                    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;
4277                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
4278                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
4279                                }
4280                            }
4281                        }
4282                        ///////////////
4283                        // process y //
4284                        ///////////////
4285                        if (add_EM_F) {
4286                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4287                            for (index_t k=0; k<numEq; k++) {
4288                                const double tmp0_1 = w3*y_p[k];
4289                                EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
4290                                EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
4291                            }
4292                        }
4293                        const index_t firstNode=m_N0*(m_N1-2)+k0;
4294                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4295                                firstNode, numEq, numComp);
4296                    }
4297                } // end colouring
4298            }
4299        } // end of parallel section
4300  }  }
4301    
4302    

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

  ViewVC Help
Powered by ViewVC 1.1.26