/[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

branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp revision 3699 by caltinay, Thu Dec 1 22:59:42 2011 UTC trunk/ripley/src/Rectangle.cpp revision 3806 by caltinay, Mon Feb 6 02:32:48 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2011 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 13  Line 13 
13    
14  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
15  extern "C" {  extern "C" {
16  #include "paso/SystemMatrixPattern.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 43  Rectangle::Rectangle(int n0, int n1, dou Line 46  Rectangle::Rectangle(int n0, int n1, dou
46      if (m_NX*m_NY != m_mpiInfo->size)      if (m_NX*m_NY != m_mpiInfo->size)
47          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
48    
49      if (n0%m_NX > 0 || n1%m_NY > 0)      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)
50          throw RipleyException("Number of elements must be separable into number of ranks in each dimension");          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");
51    
52      // local number of elements      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
53      m_NE0 = n0/m_NX;          throw RipleyException("Too few elements for the number of ranks");
54      m_NE1 = n1/m_NY;  
55      // local number of nodes (not necessarily owned)      // local number of elements (with and without overlap)
56        m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
57        if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
58            m_NE0++;
59        else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)
60            m_ownNE0--;
61    
62        m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
63        if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)
64            m_NE1++;
65        else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)
66            m_ownNE1--;
67    
68        // local number of nodes
69      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
70      m_N1 = m_NE1+1;      m_N1 = m_NE1+1;
71    
72      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
73      m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);
74      m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);      if (m_offset0 > 0)
75            m_offset0--;
76        m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);
77        if (m_offset1 > 0)
78            m_offset1--;
79    
80      populateSampleIds();      populateSampleIds();
81        createPattern();
82  }  }
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 69  string Rectangle::getDescription() const Line 94  string Rectangle::getDescription() const
94    
95  bool Rectangle::operator==(const AbstractDomain& other) const  bool Rectangle::operator==(const AbstractDomain& other) const
96  {  {
97      if (dynamic_cast<const Rectangle*>(&other))      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
98          return this==&other;      if (o) {
99            return (RipleyDomain::operator==(other) &&
100                    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
103                    && m_NX==o->m_NX && m_NY==o->m_NY);
104        }
105    
106      return false;      return false;
107  }  }
# Line 83  void Rectangle::dump(const string& fileN Line 114  void Rectangle::dump(const string& fileN
114          fn+=".silo";          fn+=".silo";
115      }      }
116    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
117      int driver=DB_HDF5;          int driver=DB_HDF5;    
     string siloPath;  
118      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
119        const char* blockDirFmt = "/block%04d";
120    
121  #ifdef ESYS_MPI  #ifdef ESYS_MPI
122      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
123        const int NUM_SILO_FILES = 1;
124  #endif  #endif
125    
126      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 106  void Rectangle::dump(const string& fileN Line 136  void Rectangle::dump(const string& fileN
136                          PMPIO_DefaultClose, (void*)&driver);                          PMPIO_DefaultClose, (void*)&driver);
137          }          }
138          if (baton) {          if (baton) {
139              char str[64];              char siloPath[64];
140              snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));              snprintf(siloPath, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
141              siloPath = str;              dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath);
             dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());  
142          }          }
143  #endif  #endif
144      } else {      } else {
# Line 121  void Rectangle::dump(const string& fileN Line 150  void Rectangle::dump(const string& fileN
150              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
151                      getDescription().c_str(), driver);                      getDescription().c_str(), driver);
152          }          }
153            char siloPath[64];
154            snprintf(siloPath, 64, blockDirFmt, 0);
155            DBMkDir(dbfile, siloPath);
156            DBSetDir(dbfile, siloPath);
157      }      }
158    
159      if (!dbfile)      if (!dbfile)
# Line 142  void Rectangle::dump(const string& fileN Line 175  void Rectangle::dump(const string& fileN
175      pair<double,double> ydy = getFirstCoordAndSpacing(1);      pair<double,double> ydy = getFirstCoordAndSpacing(1);
176  #pragma omp parallel  #pragma omp parallel
177      {      {
178  #pragma omp for  #pragma omp for nowait
179          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_N0; i0++) {
180              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=xdx.first+i0*xdx.second;
181          }          }
182  #pragma omp for  #pragma omp for nowait
183          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_N1; i1++) {
184              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=ydy.first+i1*ydy.second;
185          }          }
# Line 213  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 221  const int* Rectangle::borrowSampleRefere Line 254  const int* Rectangle::borrowSampleRefere
254  {  {
255      switch (fsType) {      switch (fsType) {
256          case Nodes:          case Nodes:
257            case ReducedNodes: // FIXME: reduced
258              return &m_nodeId[0];              return &m_nodeId[0];
259            case DegreesOfFreedom:
260            case ReducedDegreesOfFreedom: // FIXME: reduced
261                return &m_dofId[0];
262          case Elements:          case Elements:
263            case ReducedElements:
264              return &m_elementId[0];              return &m_elementId[0];
265          case FaceElements:          case FaceElements:
266            case ReducedFaceElements:
267              return &m_faceId[0];              return &m_faceId[0];
268          default:          default:
269              break;              break;
270      }      }
271    
272      stringstream msg;      stringstream msg;
273      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs: invalid function space type " << fsType;
         << fsType;  
274      throw RipleyException(msg.str());      throw RipleyException(msg.str());
275  }  }
276    
277  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsType, index_t id) const
 {  
 #ifdef ESYS_MPI  
     if (fsCode == Nodes) {  
         const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];  
         const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;  
         return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);  
     } else  
         throw RipleyException("ownSample() only implemented for Nodes");  
 #else  
     return true;  
 #endif  
 }  
   
 void Rectangle::interpolateOnDomain(escript::Data& target,  
                                     const escript::Data& in) const  
278  {  {
279      const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));      if (getMPISize()==1)
280      const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));          return true;
     if (inDomain != *this)  
         throw RipleyException("Illegal domain of interpolant");  
     if (targetDomain != *this)  
         throw RipleyException("Illegal domain of interpolation target");  
281    
282      stringstream msg;      switch (fsType) {
     msg << "interpolateOnDomain() not implemented for function space "  
         << in.getFunctionSpace().getTypeCode() << " -> "  
         << target.getFunctionSpace().getTypeCode();  
   
     switch (in.getFunctionSpace().getTypeCode()) {  
283          case Nodes:          case Nodes:
284            case ReducedNodes: // FIXME: reduced
285                return (m_dofMap[id] < getNumDOF());
286          case DegreesOfFreedom:          case DegreesOfFreedom:
287              switch (target.getFunctionSpace().getTypeCode()) {          case ReducedDegreesOfFreedom:
288                  case DegreesOfFreedom:              return true;
289                      target=in;          case Elements:
290                      break;          case ReducedElements:
291                // check ownership of element's bottom left node
292                  case Elements:              return (m_dofMap[id%m_NE0+m_N0*(id/m_NE0)] < getNumDOF());
293                  {          case FaceElements:
294                      const double tmp0_2 = 0.62200846792814621559;          case ReducedFaceElements:
295                      const double tmp0_1 = 0.044658198738520451079;              {
296                      const double tmp0_0 = 0.16666666666666666667;                  // determine which face the sample belongs to before
297                      const dim_t numComp = in.getDataPointSize();                  // checking ownership of corresponding element's first node
298                      escript::Data* inn=const_cast<escript::Data*>(&in);                  const IndexVector faces = getNumFacesPerBoundary();
299  #pragma omp parallel for                  dim_t n=0;
300                      for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (size_t i=0; i<faces.size(); i++) {
301                          for (index_t k0=0; k0 < m_NE0; ++k0) {                      n+=faces[i];
302                              const register double* f_10 = inn->getSampleDataRO(INDEX2(k0+1,k1, m_N0));                      if (id<n) {
303                              const register double* f_11 = inn->getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                          index_t k;
304                              const register double* f_01 = inn->getSampleDataRO(INDEX2(k0,k1+1, m_N0));                          if (i==1)
305                              const register double* f_00 = inn->getSampleDataRO(INDEX2(k0,k1, m_N0));                              k=m_N0-2;
306                              double* o = target.getSampleDataRW(INDEX2(k0,k1,m_NE0));                          else if (i==3)
307                              for (index_t i=0; i < numComp; ++i) {                              k=m_N0*(m_N1-2);
308                                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);                          else
309                                  o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);                              k=0;
310                                  o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);                          // determine whether to move right or up
311                                  o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);                          const index_t delta=(i/2==0 ? m_N0 : 1);
312                              } // close component loop i                          return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());
313                          } // close k0 loop                      }
                     } // close k1 loop  
                     break;  
314                  }                  }
315                    return false;
                 default:  
                     throw RipleyException(msg.str());  
316              }              }
             break;  
317          default:          default:
318              throw RipleyException(msg.str());              break;
319      }      }
320    
321        stringstream msg;
322        msg << "ownSample: invalid function space type " << fsType;
323        throw RipleyException(msg.str());
324  }  }
325    
326  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  void Rectangle::setToNormal(escript::Data& out) const
                                                 bool reducedColOrder) const  
327  {  {
328      if (reducedRowOrder || reducedColOrder)      if (out.getFunctionSpace().getTypeCode() == FaceElements) {
329          throw RipleyException("getPattern() not implemented for reduced order");          out.requireWrite();
330    #pragma omp parallel
331            {
332                if (m_faceOffset[0] > -1) {
333    #pragma omp for nowait
334                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
335                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
336                        // set vector at two quadrature points
337                        *o++ = -1.;
338                        *o++ = 0.;
339                        *o++ = -1.;
340                        *o = 0.;
341                    }
342                }
343    
344      // connector              if (m_faceOffset[1] > -1) {
345      RankVector neighbour;  #pragma omp for nowait
346      IndexVector offsetInShared(1,0);                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
347      IndexVector sendShared, recvShared;                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
348      const IndexVector faces=getNumFacesPerBoundary();                      // set vector at two quadrature points
349      const index_t left = (m_offset0==0 ? 0 : 1);                      *o++ = 1.;
350      const index_t bottom = (m_offset1==0 ? 0 : 1);                      *o++ = 0.;
351      // corner node from bottom-left                      *o++ = 1.;
352      if (faces[0] == 0 && faces[2] == 0) {                      *o = 0.;
353          neighbour.push_back(m_mpiInfo->rank-m_NX-1);                  }
354          offsetInShared.push_back(offsetInShared.back()+1);              }
         sendShared.push_back(m_nodeId[m_N0+left]);  
         recvShared.push_back(m_nodeId[0]);  
     }  
     // bottom edge  
     if (faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX);  
         offsetInShared.push_back(offsetInShared.back()+m_N0-left);  
         for (dim_t i=left; i<m_N0; i++) {  
             // easy case, we know the neighbour id's  
             sendShared.push_back(m_nodeId[m_N0+i]);  
             recvShared.push_back(m_nodeId[i]);  
         }  
     }  
     // corner node from bottom-right  
     if (faces[1] == 0 && faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX+1);  
         const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour.back()/m_NX == 0 ? m_N1 : m_N1-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[(bottom+1)*m_N0-1]);  
         recvShared.push_back(first+N0*(N1-1));  
     }  
     // left edge  
     if (faces[0] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-1);  
         offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);  
         for (dim_t i=bottom; i<m_N1; i++) {  
             // easy case, we know the neighbour id's  
             sendShared.push_back(m_nodeId[i*m_N0+1]);  
             recvShared.push_back(m_nodeId[i*m_N0]);  
         }  
     }  
     // right edge  
     if (faces[1] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+1);  
         const index_t rightN0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);  
         for (dim_t i=bottom; i<m_N1; i++) {  
             sendShared.push_back(m_nodeId[(i+1)*m_N0-1]);  
             recvShared.push_back(first+rightN0*(i-bottom));  
         }  
     }  
     // corner node from top-left  
     if (faces[0] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX-1);  
         const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+left]);  
         recvShared.push_back(first+N0-1);  
     }  
     // top edge  
     if (faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+m_N0-left);  
         for (dim_t i=left; i<m_N0; i++) {  
             sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+i]);  
             recvShared.push_back(first+i-left);  
         }  
     }  
     // corner node from top-right  
     if (faces[1] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX+1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0*m_N1-1]);  
         recvShared.push_back(first);  
     }  
     const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];  
     cout << "--- rcv_shcomp ---" << endl;  
     cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;  
     for (size_t i=0; i<neighbour.size(); i++) {  
         cout << "neighbor[" << i << "]=" << neighbour[i]  
             << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;  
     }  
     for (size_t i=0; i<recvShared.size(); i++) {  
         cout << "shared[" << i << "]=" << recvShared[i] << endl;  
     }  
     cout << "--- snd_shcomp ---" << endl;  
     for (size_t i=0; i<sendShared.size(); i++) {  
         cout << "shared[" << i << "]=" << sendShared[i] << endl;  
     }  
355    
356      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(              if (m_faceOffset[2] > -1) {
357              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],  #pragma omp for nowait
358              &offsetInShared[0], 1, 0, m_mpiInfo);                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
359      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
360              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],                      // set vector at two quadrature points
361              &offsetInShared[0], 1, 0, m_mpiInfo);                      *o++ = 0.;
362      Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);                      *o++ = -1.;
363      Paso_SharedComponents_free(snd_shcomp);                      *o++ = 0.;
364      Paso_SharedComponents_free(rcv_shcomp);                      *o = -1.;
365                    }
366                }
367    
368      // create patterns              if (m_faceOffset[3] > -1) {
369      dim_t M, N;  #pragma omp for nowait
370      IndexVector ptr(1,0);                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
371      IndexVector index;                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
372                        // set vector at two quadrature points
373      // main pattern                      *o++ = 0.;
374      for (index_t i=0; i<numDOF; i++) {                      *o++ = 1.;
375          // always add the node itself                      *o++ = 0.;
376          index.push_back(i);                      *o = 1.;
377          int num=insertNeighbours(index, i);                  }
378          ptr.push_back(ptr.back()+num+1);              }
379      }          } // end of parallel section
380      M=N=ptr.size()-1;      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
381      // paso will manage the memory          out.requireWrite();
382      index_t* indexC = MEMALLOC(index.size(),index_t);  #pragma omp parallel
383      index_t* ptrC = MEMALLOC(ptr.size(), index_t);          {
384      copy(index.begin(), index.end(), indexC);              if (m_faceOffset[0] > -1) {
385      copy(ptr.begin(), ptr.end(), ptrC);  #pragma omp for nowait
386      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
387              M, N, ptrC, indexC);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
388                        *o++ = -1.;
389                        *o = 0.;
390                    }
391                }
392    
393      cout << "--- main_pattern ---" << endl;              if (m_faceOffset[1] > -1) {
394      cout << "M=" << M << ", N=" << N << endl;  #pragma omp for nowait
395      for (size_t i=0; i<ptr.size(); i++) {                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
396          cout << "ptr[" << i << "]=" << ptr[i] << endl;                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
397      }                      *o++ = 1.;
398      for (size_t i=0; i<index.size(); i++) {                      *o = 0.;
399          cout << "index[" << i << "]=" << index[i] << endl;                  }
400                }
401    
402                if (m_faceOffset[2] > -1) {
403    #pragma omp for nowait
404                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
405                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
406                        *o++ = 0.;
407                        *o = -1.;
408                    }
409                }
410    
411                if (m_faceOffset[3] > -1) {
412    #pragma omp for nowait
413                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
414                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
415                        *o++ = 0.;
416                        *o = 1.;
417                    }
418                }
419            } // end of parallel section
420    
421        } else {
422            stringstream msg;
423            msg << "setToNormal: invalid function space type "
424                << out.getFunctionSpace().getTypeCode();
425            throw RipleyException(msg.str());
426      }      }
427    }
428    
429      ptr.clear();  void Rectangle::setToSize(escript::Data& out) const
430      index.clear();  {
431        if (out.getFunctionSpace().getTypeCode() == Elements
432                || out.getFunctionSpace().getTypeCode() == ReducedElements) {
433            out.requireWrite();
434            const dim_t numQuad=out.getNumDataPointsPerSample();
435            const double hSize=getFirstCoordAndSpacing(0).second;
436            const double vSize=getFirstCoordAndSpacing(1).second;
437            const double size=sqrt(hSize*hSize+vSize*vSize);
438    #pragma omp parallel for
439            for (index_t k = 0; k < getNumElements(); ++k) {
440                double* o = out.getSampleDataRW(k);
441                fill(o, o+numQuad, size);
442            }
443        } else if (out.getFunctionSpace().getTypeCode() == FaceElements
444                || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
445            out.requireWrite();
446            const dim_t numQuad=out.getNumDataPointsPerSample();
447            const double hSize=getFirstCoordAndSpacing(0).second;
448            const double vSize=getFirstCoordAndSpacing(1).second;
449    #pragma omp parallel
450            {
451                if (m_faceOffset[0] > -1) {
452    #pragma omp for nowait
453                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
454                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
455                        fill(o, o+numQuad, vSize);
456                    }
457                }
458    
459      // column & row couple patterns              if (m_faceOffset[1] > -1) {
460      Paso_Pattern *colCouplePattern, *rowCouplePattern;  #pragma omp for nowait
461      generateCouplePatterns(&colCouplePattern, &rowCouplePattern);                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
462                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
463                        fill(o, o+numQuad, vSize);
464                    }
465                }
466    
467      // allocate paso distribution              if (m_faceOffset[2] > -1) {
468      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,  #pragma omp for nowait
469              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
470                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
471                        fill(o, o+numQuad, hSize);
472                    }
473                }
474    
475      Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(              if (m_faceOffset[3] > -1) {
476              MATRIX_FORMAT_DEFAULT, distribution, distribution,  #pragma omp for nowait
477              mainPattern, colCouplePattern, rowCouplePattern,                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
478              connector, connector);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
479      Paso_Pattern_free(mainPattern);                      fill(o, o+numQuad, hSize);
480      Paso_Pattern_free(colCouplePattern);                  }
481      Paso_Pattern_free(rowCouplePattern);              }
482      Paso_Distribution_free(distribution);          } // end of parallel section
483      return pattern;  
484        } else {
485            stringstream msg;
486            msg << "setToSize: invalid function space type "
487                << out.getFunctionSpace().getTypeCode();
488            throw RipleyException(msg.str());
489        }
490    }
491    
492    Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
493                                                    bool reducedColOrder) const
494    {
495        /* FIXME: reduced
496        if (reducedRowOrder || reducedColOrder)
497            throw RipleyException("getPattern() not implemented for reduced order");
498        */
499        return m_pattern;
500  }  }
501    
502  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
# Line 524  IndexVector Rectangle::getNumFacesPerBou Line 550  IndexVector Rectangle::getNumFacesPerBou
550      return ret;      return ret;
551  }  }
552    
553    IndexVector Rectangle::getNumSubdivisionsPerDim() const
554    {
555        IndexVector ret;
556        ret.push_back(m_NX);
557        ret.push_back(m_NY);
558        return ret;
559    }
560    
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
572    dim_t Rectangle::getNumDOF() const
573    {
574        return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;
575  }  }
576    
577  //protected  //protected
# Line 567  void Rectangle::assembleCoordinates(escr Line 607  void Rectangle::assembleCoordinates(escr
607      }      }
608  }  }
609    
610    //protected
611    void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
612    {
613        const dim_t numComp = in.getDataPointSize();
614        const double h0 = m_l0/m_gNE0;
615        const double h1 = m_l1/m_gNE1;
616        const double cx0 = -1./h0;
617        const double cx1 = -.78867513459481288225/h0;
618        const double cx2 = -.5/h0;
619        const double cx3 = -.21132486540518711775/h0;
620        const double cx4 = .21132486540518711775/h0;
621        const double cx5 = .5/h0;
622        const double cx6 = .78867513459481288225/h0;
623        const double cx7 = 1./h0;
624        const double cy0 = -1./h1;
625        const double cy1 = -.78867513459481288225/h1;
626        const double cy2 = -.5/h1;
627        const double cy3 = -.21132486540518711775/h1;
628        const double cy4 = .21132486540518711775/h1;
629        const double cy5 = .5/h1;
630        const double cy6 = .78867513459481288225/h1;
631        const double cy7 = 1./h1;
632    
633        if (out.getFunctionSpace().getTypeCode() == Elements) {
634            out.requireWrite();
635    #pragma omp parallel for
636            for (index_t k1=0; k1 < m_NE1; ++k1) {
637                for (index_t k0=0; k0 < m_NE0; ++k0) {
638                    const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
639                    const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
640                    const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
641                    const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
642                    double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
643                    for (index_t i=0; i < numComp; ++i) {
644                        o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
645                        o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
646                        o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
647                        o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
648                        o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
649                        o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
650                        o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
651                        o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
652                    } // end of component loop i
653                } // end of k0 loop
654            } // end of k1 loop
655        } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
656            out.requireWrite();
657    #pragma omp parallel for
658            for (index_t k1=0; k1 < m_NE1; ++k1) {
659                for (index_t k0=0; k0 < m_NE0; ++k0) {
660                    const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
661                    const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
662                    const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
663                    const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
664                    double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
665                    for (index_t i=0; i < numComp; ++i) {
666                        o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
667                        o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
668                    } // end of component loop i
669                } // end of k0 loop
670            } // end of k1 loop
671    
672        } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
673            out.requireWrite();
674    #pragma omp parallel
675            {
676                if (m_faceOffset[0] > -1) {
677    #pragma omp for nowait
678                    for (index_t k1=0; k1 < m_NE1; ++k1) {
679                        const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
680                        const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
681                        const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
682                        const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
683                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
684                        for (index_t i=0; i < numComp; ++i) {
685                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
686                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
687                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
688                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
689                        } // end of component loop i
690                    } // end of k1 loop
691                } // end of face 0
692                if (m_faceOffset[1] > -1) {
693    #pragma omp for nowait
694                    for (index_t k1=0; k1 < m_NE1; ++k1) {
695                        const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
696                        const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
697                        const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
698                        const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
699                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
700                        for (index_t i=0; i < numComp; ++i) {
701                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
702                            o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
703                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
704                            o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
705                        } // end of component loop i
706                    } // end of k1 loop
707                } // end of face 1
708                if (m_faceOffset[2] > -1) {
709    #pragma omp for nowait
710                    for (index_t k0=0; k0 < m_NE0; ++k0) {
711                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
712                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
713                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
714                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
715                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
716                        for (index_t i=0; i < numComp; ++i) {
717                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
718                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
719                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
720                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
721                        } // end of component loop i
722                    } // end of k0 loop
723                } // end of face 2
724                if (m_faceOffset[3] > -1) {
725    #pragma omp for nowait
726                    for (index_t k0=0; k0 < m_NE0; ++k0) {
727                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
728                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
729                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
730                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
731                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
732                        for (index_t i=0; i < numComp; ++i) {
733                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
734                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
735                            o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
736                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
737                        } // end of component loop i
738                    } // end of k0 loop
739                } // end of face 3
740            } // end of parallel section
741    
742        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
743            out.requireWrite();
744    #pragma omp parallel
745            {
746                if (m_faceOffset[0] > -1) {
747    #pragma omp for nowait
748                    for (index_t k1=0; k1 < m_NE1; ++k1) {
749                        const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
750                        const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
751                        const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
752                        const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
753                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
754                        for (index_t i=0; i < numComp; ++i) {
755                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
756                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
757                        } // end of component loop i
758                    } // end of k1 loop
759                } // end of face 0
760                if (m_faceOffset[1] > -1) {
761    #pragma omp for nowait
762                    for (index_t k1=0; k1 < m_NE1; ++k1) {
763                        const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
764                        const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
765                        const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
766                        const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
767                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
768                        for (index_t i=0; i < numComp; ++i) {
769                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
770                            o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
771                        } // end of component loop i
772                    } // end of k1 loop
773                } // end of face 1
774                if (m_faceOffset[2] > -1) {
775    #pragma omp for nowait
776                    for (index_t k0=0; k0 < m_NE0; ++k0) {
777                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
778                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
779                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
780                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
781                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
782                        for (index_t i=0; i < numComp; ++i) {
783                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
784                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
785                        } // end of component loop i
786                    } // end of k0 loop
787                } // end of face 2
788                if (m_faceOffset[3] > -1) {
789    #pragma omp for nowait
790                    for (index_t k0=0; k0 < m_NE0; ++k0) {
791                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
792                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
793                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
794                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
795                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
796                        for (index_t i=0; i < numComp; ++i) {
797                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
798                            o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
799                        } // end of component loop i
800                    } // end of k0 loop
801                } // end of face 3
802            } // end of parallel section
803        }
804    }
805    
806    //protected
807    void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
808    {
809        const dim_t numComp = arg.getDataPointSize();
810        const double h0 = m_l0/m_gNE0;
811        const double h1 = m_l1/m_gNE1;
812        const index_t left = (m_offset0==0 ? 0 : 1);
813        const index_t bottom = (m_offset1==0 ? 0 : 1);
814        const int fs=arg.getFunctionSpace().getTypeCode();
815        if (fs == Elements && arg.actsExpanded()) {
816    #pragma omp parallel
817            {
818                vector<double> int_local(numComp, 0);
819                const double w = h0*h1/4.;
820    #pragma omp for nowait
821                for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
822                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
823                        const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
824                        for (index_t i=0; i < numComp; ++i) {
825                            const double f0 = f[INDEX2(i,0,numComp)];
826                            const double f1 = f[INDEX2(i,1,numComp)];
827                            const double f2 = f[INDEX2(i,2,numComp)];
828                            const double f3 = f[INDEX2(i,3,numComp)];
829                            int_local[i]+=(f0+f1+f2+f3)*w;
830                        }  // end of component loop i
831                    } // end of k0 loop
832                } // end of k1 loop
833    #pragma omp critical
834                for (index_t i=0; i<numComp; i++)
835                    integrals[i]+=int_local[i];
836            } // end of parallel section
837    
838        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
839            const double w = h0*h1;
840    #pragma omp parallel
841            {
842                vector<double> int_local(numComp, 0);
843    #pragma omp for nowait
844                for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
845                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
846                        const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
847                        for (index_t i=0; i < numComp; ++i) {
848                            int_local[i]+=f[i]*w;
849                        }
850                    }
851                }
852    #pragma omp critical
853                for (index_t i=0; i<numComp; i++)
854                    integrals[i]+=int_local[i];
855            } // end of parallel section
856    
857        } else if (fs == FaceElements && arg.actsExpanded()) {
858    #pragma omp parallel
859            {
860                vector<double> int_local(numComp, 0);
861                const double w0 = h0/2.;
862                const double w1 = h1/2.;
863                if (m_faceOffset[0] > -1) {
864    #pragma omp for nowait
865                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
866                        const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
867                        for (index_t i=0; i < numComp; ++i) {
868                            const double f0 = f[INDEX2(i,0,numComp)];
869                            const double f1 = f[INDEX2(i,1,numComp)];
870                            int_local[i]+=(f0+f1)*w1;
871                        }  // end of component loop i
872                    } // end of k1 loop
873                }
874    
875                if (m_faceOffset[1] > -1) {
876    #pragma omp for nowait
877                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
878                        const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
879                        for (index_t i=0; i < numComp; ++i) {
880                            const double f0 = f[INDEX2(i,0,numComp)];
881                            const double f1 = f[INDEX2(i,1,numComp)];
882                            int_local[i]+=(f0+f1)*w1;
883                        }  // end of component loop i
884                    } // end of k1 loop
885                }
886    
887                if (m_faceOffset[2] > -1) {
888    #pragma omp for nowait
889                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
890                        const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
891                        for (index_t i=0; i < numComp; ++i) {
892                            const double f0 = f[INDEX2(i,0,numComp)];
893                            const double f1 = f[INDEX2(i,1,numComp)];
894                            int_local[i]+=(f0+f1)*w0;
895                        }  // end of component loop i
896                    } // end of k0 loop
897                }
898    
899                if (m_faceOffset[3] > -1) {
900    #pragma omp for nowait
901                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
902                        const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
903                        for (index_t i=0; i < numComp; ++i) {
904                            const double f0 = f[INDEX2(i,0,numComp)];
905                            const double f1 = f[INDEX2(i,1,numComp)];
906                            int_local[i]+=(f0+f1)*w0;
907                        }  // end of component loop i
908                    } // end of k0 loop
909                }
910    #pragma omp critical
911                for (index_t i=0; i<numComp; i++)
912                    integrals[i]+=int_local[i];
913            } // end of parallel section
914    
915        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
916    #pragma omp parallel
917            {
918                vector<double> int_local(numComp, 0);
919                if (m_faceOffset[0] > -1) {
920    #pragma omp for nowait
921                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
922                        const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
923                        for (index_t i=0; i < numComp; ++i) {
924                            int_local[i]+=f[i]*h1;
925                        }
926                    }
927                }
928    
929                if (m_faceOffset[1] > -1) {
930    #pragma omp for nowait
931                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
932                        const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
933                        for (index_t i=0; i < numComp; ++i) {
934                            int_local[i]+=f[i]*h1;
935                        }
936                    }
937                }
938    
939                if (m_faceOffset[2] > -1) {
940    #pragma omp for nowait
941                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
942                        const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
943                        for (index_t i=0; i < numComp; ++i) {
944                            int_local[i]+=f[i]*h0;
945                        }
946                    }
947                }
948    
949                if (m_faceOffset[3] > -1) {
950    #pragma omp for nowait
951                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
952                        const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
953                        for (index_t i=0; i < numComp; ++i) {
954                            int_local[i]+=f[i]*h0;
955                        }
956                    }
957                }
958    
959    #pragma omp critical
960                for (index_t i=0; i<numComp; i++)
961                    integrals[i]+=int_local[i];
962            } // end of parallel section
963        } // function space selector
964    }
965    
966    //protected
967    dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
968    {
969        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
970        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
971        const int x=node%nDOF0;
972        const int y=node/nDOF0;
973        dim_t num=0;
974        // loop through potential neighbours and add to index if positions are
975        // within bounds
976        for (int i1=-1; i1<2; i1++) {
977            for (int i0=-1; i0<2; i0++) {
978                // skip node itself
979                if (i0==0 && i1==0)
980                    continue;
981                // location of neighbour node
982                const int nx=x+i0;
983                const int ny=y+i1;
984                if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {
985                    index.push_back(ny*nDOF0+nx);
986                    num++;
987                }
988            }
989        }
990    
991        return num;
992    }
993    
994    //protected
995    void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const
996    {
997        const dim_t numComp = in.getDataPointSize();
998        out.requireWrite();
999    
1000        const index_t left = (m_offset0==0 ? 0 : 1);
1001        const index_t bottom = (m_offset1==0 ? 0 : 1);
1002        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1003        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1004    #pragma omp parallel for
1005        for (index_t i=0; i<nDOF1; i++) {
1006            for (index_t j=0; j<nDOF0; j++) {
1007                const index_t n=j+left+(i+bottom)*m_N0;
1008                const double* src=in.getSampleDataRO(n);
1009                copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
1010            }
1011        }
1012    }
1013    
1014    //protected
1015    void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const
1016    {
1017        const dim_t numComp = in.getDataPointSize();
1018        Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
1019        in.requireWrite();
1020        Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));
1021    
1022        const dim_t numDOF = getNumDOF();
1023        out.requireWrite();
1024        const double* buffer = Paso_Coupler_finishCollect(coupler);
1025    
1026    #pragma omp parallel for
1027        for (index_t i=0; i<getNumNodes(); i++) {
1028            const double* src=(m_dofMap[i]<numDOF ?
1029                    in.getSampleDataRO(m_dofMap[i])
1030                    : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1031            copy(src, src+numComp, out.getSampleDataRW(i));
1032        }
1033    }
1034    
1035  //private  //private
1036  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
1037  {  {
1038      // identifiers are ordered from left to right, bottom to top on each rank,      // identifiers are ordered from left to right, bottom to top globablly.
     // except for the shared nodes which are owned by the rank below / to the  
     // left of the current rank  
1039    
1040      // build node distribution vector first.      // build node distribution vector first.
     // m_nodeDistribution[i] is the first node id on rank i, that is  
1041      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
1042      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1043      m_nodeDistribution[1]=getNumNodes();      const dim_t numDOF=getNumDOF();
1044      for (dim_t k=1; k<m_mpiInfo->size-1; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
1045          const index_t x=k%m_NX;          m_nodeDistribution[k]=k*numDOF;
         const index_t y=k/m_NX;  
         index_t numNodes=getNumNodes();  
         if (x>0)  
             numNodes-=m_N1;  
         if (y>0)  
             numNodes-=m_N0;  
         if (x>0 && y>0)  
             numNodes++; // subtracted corner twice -> fix that  
         m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;  
1046      }      }
1047      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
   
1048      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
1049        m_dofId.resize(numDOF);
1050        m_elementId.resize(getNumElements());
1051        m_faceId.resize(getNumFaceElements());
1052    
1053    #pragma omp parallel
1054        {
1055            // nodes
1056    #pragma omp for nowait
1057            for (dim_t i1=0; i1<m_N1; i1++) {
1058                for (dim_t i0=0; i0<m_N0; i0++) {
1059                    m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;
1060                }
1061            }
1062    
1063      // the bottom row and left column are not owned by this rank so the          // degrees of freedom
1064      // identifiers need to be computed accordingly  #pragma omp for nowait
1065            for (dim_t k=0; k<numDOF; k++)
1066                m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;
1067    
1068            // elements
1069    #pragma omp for nowait
1070            for (dim_t i1=0; i1<m_NE1; i1++) {
1071                for (dim_t i0=0; i0<m_NE0; i0++) {
1072                    m_elementId[i0+i1*m_NE0]=(m_offset1+i1)*m_gNE0+m_offset0+i0;
1073                }
1074            }
1075    
1076            // face elements
1077    #pragma omp for
1078            for (dim_t k=0; k<getNumFaceElements(); k++)
1079                m_faceId[k]=k;
1080        } // end parallel section
1081    
1082        m_nodeTags.assign(getNumNodes(), 0);
1083        updateTagsInUse(Nodes);
1084    
1085        m_elementTags.assign(getNumElements(), 0);
1086        updateTagsInUse(Elements);
1087    
1088        // generate face offset vector and set face tags
1089        const IndexVector facesPerEdge = getNumFacesPerBoundary();
1090        const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
1091        const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
1092        m_faceOffset.assign(facesPerEdge.size(), -1);
1093        m_faceTags.clear();
1094        index_t offset=0;
1095        for (size_t i=0; i<facesPerEdge.size(); i++) {
1096            if (facesPerEdge[i]>0) {
1097                m_faceOffset[i]=offset;
1098                offset+=facesPerEdge[i];
1099                m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);
1100            }
1101        }
1102        setTagMap("left", LEFT);
1103        setTagMap("right", RIGHT);
1104        setTagMap("bottom", BOTTOM);
1105        setTagMap("top", TOP);
1106        updateTagsInUse(FaceElements);
1107    }
1108    
1109    //private
1110    void Rectangle::createPattern()
1111    {
1112        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1113        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1114      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1115      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1116      if (left>0) {  
1117          const int neighbour=m_mpiInfo->rank-1;      // populate node->DOF mapping with own degrees of freedom.
1118          const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      // The rest is assigned in the loop further down
1119        m_dofMap.assign(getNumNodes(), 0);
1120  #pragma omp parallel for  #pragma omp parallel for
1121          for (dim_t i1=bottom; i1<m_N1; i1++) {      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1122              m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]          for (index_t j=left; j<left+nDOF0; j++) {
1123                  + (i1-bottom+1)*leftN0-1;              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
1124          }          }
1125      }      }
1126      if (bottom>0) {  
1127          const int neighbour=m_mpiInfo->rank-m_NX;      // build list of shared components and neighbours by looping through
1128          const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      // all potential neighbouring ranks and checking if positions are
1129          const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);      // within bounds
1130  #pragma omp parallel for      const dim_t numDOF=nDOF0*nDOF1;
1131          for (dim_t i0=left; i0<m_N0; i0++) {      vector<IndexVector> colIndices(numDOF); // for the couple blocks
1132              m_nodeId[i0]=m_nodeDistribution[neighbour]      RankVector neighbour;
1133                  + (bottomN1-1)*bottomN0 + i0 - left;      IndexVector offsetInShared(1,0);
1134        IndexVector sendShared, recvShared;
1135        int numShared=0;
1136        const int x=m_mpiInfo->rank%m_NX;
1137        const int y=m_mpiInfo->rank/m_NX;
1138        for (int i1=-1; i1<2; i1++) {
1139            for (int i0=-1; i0<2; i0++) {
1140                // skip this rank
1141                if (i0==0 && i1==0)
1142                    continue;
1143                // location of neighbour rank
1144                const int nx=x+i0;
1145                const int ny=y+i1;
1146                if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {
1147                    neighbour.push_back(ny*m_NX+nx);
1148                    if (i0==0) {
1149                        // sharing top or bottom edge
1150                        const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1151                        const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);
1152                        offsetInShared.push_back(offsetInShared.back()+nDOF0);
1153                        for (dim_t i=0; i<nDOF0; i++, numShared++) {
1154                            sendShared.push_back(firstDOF+i);
1155                            recvShared.push_back(numDOF+numShared);
1156                            if (i>0)
1157                                colIndices[firstDOF+i-1].push_back(numShared);
1158                            colIndices[firstDOF+i].push_back(numShared);
1159                            if (i<nDOF0-1)
1160                                colIndices[firstDOF+i+1].push_back(numShared);
1161                            m_dofMap[firstNode+i]=numDOF+numShared;
1162                        }
1163                    } else if (i1==0) {
1164                        // sharing left or right edge
1165                        const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1166                        const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);
1167                        offsetInShared.push_back(offsetInShared.back()+nDOF1);
1168                        for (dim_t i=0; i<nDOF1; i++, numShared++) {
1169                            sendShared.push_back(firstDOF+i*nDOF0);
1170                            recvShared.push_back(numDOF+numShared);
1171                            if (i>0)
1172                                colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);
1173                            colIndices[firstDOF+i*nDOF0].push_back(numShared);
1174                            if (i<nDOF1-1)
1175                                colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1176                            m_dofMap[firstNode+i*m_N0]=numDOF+numShared;
1177                        }
1178                    } else {
1179                        // sharing a node
1180                        const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1181                        const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);
1182                        offsetInShared.push_back(offsetInShared.back()+1);
1183                        sendShared.push_back(dof);
1184                        recvShared.push_back(numDOF+numShared);
1185                        colIndices[dof].push_back(numShared);
1186                        m_dofMap[node]=numDOF+numShared;
1187                        ++numShared;
1188                    }
1189                }
1190          }          }
1191      }      }
1192      if (left>0 && bottom>0) {  
1193          const int neighbour=m_mpiInfo->rank-m_NX-1;      // create connector
1194          const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1195          const index_t N1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
1196          m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1-1;              &offsetInShared[0], 1, 0, m_mpiInfo);
1197        Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
1198                numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
1199                &offsetInShared[0], 1, 0, m_mpiInfo);
1200        m_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1201        Paso_SharedComponents_free(snd_shcomp);
1202        Paso_SharedComponents_free(rcv_shcomp);
1203    
1204        // create main and couple blocks
1205        Paso_Pattern *mainPattern = createMainPattern();
1206        Paso_Pattern *colPattern, *rowPattern;
1207        createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
1208    
1209        // allocate paso distribution
1210        Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
1211                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
1212    
1213        // finally create the system matrix
1214        m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
1215                distribution, distribution, mainPattern, colPattern, rowPattern,
1216                m_connector, m_connector);
1217    
1218        Paso_Distribution_free(distribution);
1219    
1220        // useful debug output
1221        /*
1222        cout << "--- rcv_shcomp ---" << endl;
1223        cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
1224        for (size_t i=0; i<neighbour.size(); i++) {
1225            cout << "neighbor[" << i << "]=" << neighbour[i]
1226                << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
1227        }
1228        for (size_t i=0; i<recvShared.size(); i++) {
1229            cout << "shared[" << i << "]=" << recvShared[i] << endl;
1230        }
1231        cout << "--- snd_shcomp ---" << endl;
1232        for (size_t i=0; i<sendShared.size(); i++) {
1233            cout << "shared[" << i << "]=" << sendShared[i] << endl;
1234        }
1235        cout << "--- dofMap ---" << endl;
1236        for (size_t i=0; i<m_dofMap.size(); i++) {
1237            cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
1238        }
1239        cout << "--- colIndices ---" << endl;
1240        for (size_t i=0; i<colIndices.size(); i++) {
1241            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
1242      }      }
1243        */
1244    
1245      // the rest of the id's are contiguous      /*
1246      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];      cout << "--- main_pattern ---" << endl;
1247  #pragma omp parallel for      cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
1248      for (dim_t i1=bottom; i1<m_N1; i1++) {      for (size_t i=0; i<mainPattern->numOutput+1; i++) {
1249          for (dim_t i0=left; i0<m_N0; i0++) {          cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
             m_nodeId[i0+i1*m_N0] = firstId+i0-left+(i1-bottom)*(m_N0-left);  
         }  
1250      }      }
1251        for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
1252            cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
1253        }
1254        */
1255    
1256      // elements      /*
1257      m_elementId.resize(getNumElements());      cout << "--- colCouple_pattern ---" << endl;
1258  #pragma omp parallel for      cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
1259      for (dim_t k=0; k<getNumElements(); k++) {      for (size_t i=0; i<colPattern->numOutput+1; i++) {
1260          m_elementId[k]=k;          cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
1261        }
1262        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
1263            cout << "index[" << i << "]=" << colPattern->index[i] << endl;
1264      }      }
1265        */
1266    
1267      // face elements      /*
1268      m_faceId.resize(getNumFaceElements());      cout << "--- rowCouple_pattern ---" << endl;
1269  #pragma omp parallel for      cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
1270      for (dim_t k=0; k<getNumFaceElements(); k++) {      for (size_t i=0; i<rowPattern->numOutput+1; i++) {
1271          m_faceId[k]=k;          cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
1272        }
1273        for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
1274            cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
1275      }      }
1276        */
1277    
1278        Paso_Pattern_free(mainPattern);
1279        Paso_Pattern_free(colPattern);
1280        Paso_Pattern_free(rowPattern);
1281  }  }
1282    
1283  //private  //private
1284  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const  void Rectangle::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
1285             const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
1286             bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
1287  {  {
1288      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);      IndexVector rowIndex;
1289      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);      rowIndex.push_back(m_dofMap[firstNode]);
1290      const int x=node%myN0;      rowIndex.push_back(m_dofMap[firstNode+1]);
1291      const int y=node/myN0;      rowIndex.push_back(m_dofMap[firstNode+m_N0]);
1292      int num=0;      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
1293      if (y>0) {      if (addF) {
1294          if (x>0) {          double *F_p=F.getSampleDataRW(0);
1295              // bottom-left          for (index_t i=0; i<rowIndex.size(); i++) {
1296              index.push_back(node-myN0-1);              if (rowIndex[i]<getNumDOF()) {
1297              num++;                  for (index_t eq=0; eq<nEq; eq++) {
1298          }                      F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
1299          // bottom                  }
1300          index.push_back(node-myN0);              }
1301          num++;          }
1302          if (x<myN0-1) {      }
1303              // bottom-right      if (addS) {
1304              index.push_back(node-myN0+1);          addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
             num++;  
         }  
     }  
     if (x<myN0-1) {  
         // right  
         index.push_back(node+1);  
         num++;  
         if (y<myN1-1) {  
             // top-right  
             index.push_back(node+myN0+1);  
             num++;  
         }  
     }  
     if (y<myN1-1) {  
         // top  
         index.push_back(node+myN0);  
         num++;  
         if (x>0) {  
             // top-left  
             index.push_back(node+myN0-1);  
             num++;  
         }  
     }  
     if (x>0) {  
         // left  
         index.push_back(node-1);  
         num++;  
1305      }      }
1306    }
1307    
1308      return num;  //protected
1309    void Rectangle::interpolateNodesOnElements(escript::Data& out,
1310                                            escript::Data& in, bool reduced) const
1311    {
1312        const dim_t numComp = in.getDataPointSize();
1313        if (reduced) {
1314            out.requireWrite();
1315            const double c0 = .25;
1316    #pragma omp parallel for
1317            for (index_t k1=0; k1 < m_NE1; ++k1) {
1318                for (index_t k0=0; k0 < m_NE0; ++k0) {
1319                    const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1320                    const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1321                    const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1322                    const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1323                    double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1324                    for (index_t i=0; i < numComp; ++i) {
1325                        o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1326                    } /* end of component loop i */
1327                } /* end of k0 loop */
1328            } /* end of k1 loop */
1329        } else {
1330            out.requireWrite();
1331            const double c0 = .16666666666666666667;
1332            const double c1 = .044658198738520451079;
1333            const double c2 = .62200846792814621559;
1334    #pragma omp parallel for
1335            for (index_t k1=0; k1 < m_NE1; ++k1) {
1336                for (index_t k0=0; k0 < m_NE0; ++k0) {
1337                    const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1338                    const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1339                    const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1340                    const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1341                    double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1342                    for (index_t i=0; i < numComp; ++i) {
1343                        o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);
1344                        o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);
1345                        o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);
1346                        o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);
1347                    } /* end of component loop i */
1348                } /* end of k0 loop */
1349            } /* end of k1 loop */
1350        }
1351  }  }
1352    
1353  //private  //protected
1354  void Rectangle::generateCouplePatterns(Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const  void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,
1355                                            bool reduced) const
1356  {  {
1357      IndexVector ptr(1,0);      const dim_t numComp = in.getDataPointSize();
1358      IndexVector index;      if (reduced) {
1359      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);          out.requireWrite();
1360      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);          const double c0 = .5;
1361      const IndexVector faces=getNumFacesPerBoundary();  #pragma omp parallel
1362            {
1363                if (m_faceOffset[0] > -1) {
1364    #pragma omp for nowait
1365                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1366                        const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1367                        const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1368                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1369                        for (index_t i=0; i < numComp; ++i) {
1370                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
1371                        } /* end of component loop i */
1372                    } /* end of k1 loop */
1373                } /* end of face 0 */
1374                if (m_faceOffset[1] > -1) {
1375    #pragma omp for nowait
1376                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1377                        const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1378                        const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1379                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1380                        for (index_t i=0; i < numComp; ++i) {
1381                            o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
1382                        } /* end of component loop i */
1383                    } /* end of k1 loop */
1384                } /* end of face 1 */
1385                if (m_faceOffset[2] > -1) {
1386    #pragma omp for nowait
1387                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1388                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1389                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1390                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1391                        for (index_t i=0; i < numComp; ++i) {
1392                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
1393                        } /* end of component loop i */
1394                    } /* end of k0 loop */
1395                } /* end of face 2 */
1396                if (m_faceOffset[3] > -1) {
1397    #pragma omp for nowait
1398                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1399                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1400                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1401                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1402                        for (index_t i=0; i < numComp; ++i) {
1403                            o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1404                        } /* end of component loop i */
1405                    } /* end of k0 loop */
1406                } /* end of face 3 */
1407            } // end of parallel section
1408        } else {
1409            out.requireWrite();
1410            const double c0 = 0.21132486540518711775;
1411            const double c1 = 0.78867513459481288225;
1412    #pragma omp parallel
1413            {
1414                if (m_faceOffset[0] > -1) {
1415    #pragma omp for nowait
1416                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1417                        const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1418                        const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1419                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1420                        for (index_t i=0; i < numComp; ++i) {
1421                            o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;
1422                            o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_01[i]*c1;
1423                        } /* end of component loop i */
1424                    } /* end of k1 loop */
1425                } /* end of face 0 */
1426                if (m_faceOffset[1] > -1) {
1427    #pragma omp for nowait
1428                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1429                        const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1430                        const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1431                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1432                        for (index_t i=0; i < numComp; ++i) {
1433                            o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;
1434                            o[INDEX2(i,numComp,1)] = f_10[i]*c0 + f_11[i]*c1;
1435                        } /* end of component loop i */
1436                    } /* end of k1 loop */
1437                } /* end of face 1 */
1438                if (m_faceOffset[2] > -1) {
1439    #pragma omp for nowait
1440                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1441                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1442                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1443                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1444                        for (index_t i=0; i < numComp; ++i) {
1445                            o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;
1446                            o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_10[i]*c1;
1447                        } /* end of component loop i */
1448                    } /* end of k0 loop */
1449                } /* end of face 2 */
1450                if (m_faceOffset[3] > -1) {
1451    #pragma omp for nowait
1452                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1453                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1454                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1455                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1456                        for (index_t i=0; i < numComp; ++i) {
1457                            o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;
1458                            o[INDEX2(i,numComp,1)] = f_01[i]*c0 + f_11[i]*c1;
1459                        } /* end of component loop i */
1460                    } /* end of k0 loop */
1461                } /* end of face 3 */
1462            } // end of parallel section
1463        }
1464    }
1465    
1466      // bottom edge  //protected
1467      dim_t n=0;  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1468      if (faces[0] == 0) {          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1469          index.push_back(2*(myN0+myN1+1));          const escript::Data& C, const escript::Data& D,
1470          index.push_back(2*(myN0+myN1+1)+1);          const escript::Data& X, const escript::Data& Y) const
1471          n+=2;  {
1472          if (faces[2] == 0) {      const double h0 = m_l0/m_gNE0;
1473              index.push_back(0);      const double h1 = m_l1/m_gNE1;
1474              index.push_back(1);      const double w0 = -0.1555021169820365539*h1/h0;
1475              index.push_back(2);      const double w1 = 0.041666666666666666667;
1476              n+=3;      const double w2 = -0.15550211698203655390;
1477          }      const double w3 = 0.041666666666666666667*h0/h1;
1478      } else if (faces[2] == 0) {      const double w4 = 0.15550211698203655390;
1479          index.push_back(1);      const double w5 = -0.041666666666666666667;
1480          index.push_back(2);      const double w6 = -0.01116454968463011277*h1/h0;
1481          n+=2;      const double w7 = 0.011164549684630112770;
1482      }      const double w8 = -0.011164549684630112770;
1483      // n=neighbours of bottom-left corner node      const double w9 = -0.041666666666666666667*h1/h0;
1484      ptr.push_back(ptr.back()+n);      const double w10 = -0.041666666666666666667*h0/h1;
1485      n=0;      const double w11 = 0.1555021169820365539*h1/h0;
1486      if (faces[2] == 0) {      const double w12 = 0.1555021169820365539*h0/h1;
1487          for (dim_t i=1; i<myN0-1; i++) {      const double w13 = 0.01116454968463011277*h0/h1;
1488              index.push_back(i);      const double w14 = 0.01116454968463011277*h1/h0;
1489              index.push_back(i+1);      const double w15 = 0.041666666666666666667*h1/h0;
1490              index.push_back(i+2);      const double w16 = -0.01116454968463011277*h0/h1;
1491              ptr.push_back(ptr.back()+3);      const double w17 = -0.1555021169820365539*h0/h1;
1492          }      const double w18 = -0.33333333333333333333*h1/h0;
1493          index.push_back(myN0-1);      const double w19 = 0.25;
1494          index.push_back(myN0);      const double w20 = -0.25;
1495          n+=2;      const double w21 = 0.16666666666666666667*h0/h1;
1496          if (faces[1] == 0) {      const double w22 = -0.16666666666666666667*h1/h0;
1497              index.push_back(myN0+1);      const double w23 = -0.16666666666666666667*h0/h1;
1498              index.push_back(myN0+2);      const double w24 = 0.33333333333333333333*h1/h0;
1499              index.push_back(myN0+3);      const double w25 = 0.33333333333333333333*h0/h1;
1500              n+=3;      const double w26 = 0.16666666666666666667*h1/h0;
1501        const double w27 = -0.33333333333333333333*h0/h1;
1502        const double w28 = -0.032861463941450536761*h1;
1503        const double w29 = -0.032861463941450536761*h0;
1504        const double w30 = -0.12264065304058601714*h1;
1505        const double w31 = -0.0023593469594139828636*h1;
1506        const double w32 = -0.008805202725216129906*h0;
1507        const double w33 = -0.008805202725216129906*h1;
1508        const double w34 = 0.032861463941450536761*h1;
1509        const double w35 = 0.008805202725216129906*h1;
1510        const double w36 = 0.008805202725216129906*h0;
1511        const double w37 = 0.0023593469594139828636*h1;
1512        const double w38 = 0.12264065304058601714*h1;
1513        const double w39 = 0.032861463941450536761*h0;
1514        const double w40 = -0.12264065304058601714*h0;
1515        const double w41 = -0.0023593469594139828636*h0;
1516        const double w42 = 0.0023593469594139828636*h0;
1517        const double w43 = 0.12264065304058601714*h0;
1518        const double w44 = -0.16666666666666666667*h1;
1519        const double w45 = -0.083333333333333333333*h0;
1520        const double w46 = 0.083333333333333333333*h1;
1521        const double w47 = 0.16666666666666666667*h1;
1522        const double w48 = 0.083333333333333333333*h0;
1523        const double w49 = -0.16666666666666666667*h0;
1524        const double w50 = 0.16666666666666666667*h0;
1525        const double w51 = -0.083333333333333333333*h1;
1526        const double w52 = 0.025917019497006092316*h0*h1;
1527        const double w53 = 0.0018607582807716854616*h0*h1;
1528        const double w54 = 0.0069444444444444444444*h0*h1;
1529        const double w55 = 0.09672363354357992482*h0*h1;
1530        const double w56 = 0.00049858867864229740201*h0*h1;
1531        const double w57 = 0.055555555555555555556*h0*h1;
1532        const double w58 = 0.027777777777777777778*h0*h1;
1533        const double w59 = 0.11111111111111111111*h0*h1;
1534        const double w60 = -0.19716878364870322056*h1;
1535        const double w61 = -0.19716878364870322056*h0;
1536        const double w62 = -0.052831216351296779436*h0;
1537        const double w63 = -0.052831216351296779436*h1;
1538        const double w64 = 0.19716878364870322056*h1;
1539        const double w65 = 0.052831216351296779436*h1;
1540        const double w66 = 0.19716878364870322056*h0;
1541        const double w67 = 0.052831216351296779436*h0;
1542        const double w68 = -0.5*h1;
1543        const double w69 = -0.5*h0;
1544        const double w70 = 0.5*h1;
1545        const double w71 = 0.5*h0;
1546        const double w72 = 0.1555021169820365539*h0*h1;
1547        const double w73 = 0.041666666666666666667*h0*h1;
1548        const double w74 = 0.01116454968463011277*h0*h1;
1549        const double w75 = 0.25*h0*h1;
1550    
1551        rhs.requireWrite();
1552    #pragma omp parallel
1553        {
1554            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1555    #pragma omp for
1556                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
1557                    for (index_t k0=0; k0<m_NE0; ++k0)  {
1558                        bool add_EM_S=false;
1559                        bool add_EM_F=false;
1560                        vector<double> EM_S(4*4, 0);
1561                        vector<double> EM_F(4, 0);
1562                        const index_t e = k0 + m_NE0*k1;
1563                        ///////////////
1564                        // process A //
1565                        ///////////////
1566                        if (!A.isEmpty()) {
1567                            add_EM_S=true;
1568                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1569                            if (A.actsExpanded()) {
1570                                const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1571                                const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1572                                const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1573                                const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1574                                const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1575                                const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1576                                const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1577                                const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1578                                const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1579                                const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1580                                const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1581                                const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1582                                const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1583                                const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1584                                const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1585                                const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1586                                const double tmp0_0 = A_01_0 + A_01_3;
1587                                const double tmp1_0 = A_00_0 + A_00_1;
1588                                const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1589                                const double tmp3_0 = A_00_2 + A_00_3;
1590                                const double tmp4_0 = A_10_1 + A_10_2;
1591                                const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1592                                const double tmp6_0 = A_01_3 + A_10_0;
1593                                const double tmp7_0 = A_01_0 + A_10_3;
1594                                const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1595                                const double tmp9_0 = A_01_0 + A_10_0;
1596                                const double tmp12_0 = A_11_0 + A_11_2;
1597                                const double tmp10_0 = A_01_3 + A_10_3;
1598                                const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
1599                                const double tmp13_0 = A_01_2 + A_10_1;
1600                                const double tmp11_0 = A_11_1 + A_11_3;
1601                                const double tmp18_0 = A_01_1 + A_10_1;
1602                                const double tmp15_0 = A_01_1 + A_10_2;
1603                                const double tmp16_0 = A_10_0 + A_10_3;
1604                                const double tmp17_0 = A_01_1 + A_01_2;
1605                                const double tmp19_0 = A_01_2 + A_10_2;
1606                                const double tmp0_1 = A_10_3*w8;
1607                                const double tmp1_1 = tmp0_0*w1;
1608                                const double tmp2_1 = A_01_1*w4;
1609                                const double tmp3_1 = tmp1_0*w0;
1610                                const double tmp4_1 = A_01_2*w7;
1611                                const double tmp5_1 = tmp2_0*w3;
1612                                const double tmp6_1 = tmp3_0*w6;
1613                                const double tmp7_1 = A_10_0*w2;
1614                                const double tmp8_1 = tmp4_0*w5;
1615                                const double tmp9_1 = tmp2_0*w10;
1616                                const double tmp14_1 = A_10_0*w8;
1617                                const double tmp23_1 = tmp3_0*w14;
1618                                const double tmp35_1 = A_01_0*w8;
1619                                const double tmp54_1 = tmp13_0*w8;
1620                                const double tmp20_1 = tmp9_0*w4;
1621                                const double tmp25_1 = tmp12_0*w12;
1622                                const double tmp44_1 = tmp7_0*w7;
1623                                const double tmp26_1 = tmp10_0*w4;
1624                                const double tmp52_1 = tmp18_0*w8;
1625                                const double tmp48_1 = A_10_1*w7;
1626                                const double tmp46_1 = A_01_3*w8;
1627                                const double tmp50_1 = A_01_0*w2;
1628                                const double tmp56_1 = tmp19_0*w8;
1629                                const double tmp19_1 = A_10_3*w2;
1630                                const double tmp47_1 = A_10_2*w4;
1631                                const double tmp16_1 = tmp3_0*w0;
1632                                const double tmp18_1 = tmp1_0*w6;
1633                                const double tmp31_1 = tmp11_0*w12;
1634                                const double tmp55_1 = tmp15_0*w2;
1635                                const double tmp39_1 = A_10_2*w7;
1636                                const double tmp11_1 = tmp6_0*w7;
1637                                const double tmp40_1 = tmp11_0*w17;
1638                                const double tmp34_1 = tmp15_0*w8;
1639                                const double tmp33_1 = tmp14_0*w5;
1640                                const double tmp24_1 = tmp11_0*w13;
1641                                const double tmp43_1 = tmp17_0*w5;
1642                                const double tmp15_1 = A_01_2*w4;
1643                                const double tmp53_1 = tmp19_0*w2;
1644                                const double tmp27_1 = tmp3_0*w11;
1645                                const double tmp32_1 = tmp13_0*w2;
1646                                const double tmp10_1 = tmp5_0*w9;
1647                                const double tmp37_1 = A_10_1*w4;
1648                                const double tmp38_1 = tmp5_0*w15;
1649                                const double tmp17_1 = A_01_1*w7;
1650                                const double tmp12_1 = tmp7_0*w4;
1651                                const double tmp22_1 = tmp10_0*w7;
1652                                const double tmp57_1 = tmp18_0*w2;
1653                                const double tmp28_1 = tmp9_0*w7;
1654                                const double tmp29_1 = tmp1_0*w14;
1655                                const double tmp51_1 = tmp11_0*w16;
1656                                const double tmp42_1 = tmp12_0*w16;
1657                                const double tmp49_1 = tmp12_0*w17;
1658                                const double tmp21_1 = tmp1_0*w11;
1659                                const double tmp45_1 = tmp6_0*w4;
1660                                const double tmp13_1 = tmp8_0*w1;
1661                                const double tmp36_1 = tmp16_0*w1;
1662                                const double tmp41_1 = A_01_3*w2;
1663                                const double tmp30_1 = tmp12_0*w13;
1664                                EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1665                                EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
1666                                EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
1667                                EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
1668                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
1669                                EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
1670                                EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
1671                                EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
1672                                EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1673                                EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
1674                                EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
1675                                EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
1676                                EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
1677                                EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
1678                                EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
1679                                EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1680                            } else { // constant data
1681                                const double A_00 = A_p[INDEX2(0,0,2)];
1682                                const double A_10 = A_p[INDEX2(1,0,2)];
1683                                const double A_01 = A_p[INDEX2(0,1,2)];
1684                                const double A_11 = A_p[INDEX2(1,1,2)];
1685                                const double tmp0_0 = A_01 + A_10;
1686                                const double tmp0_1 = A_00*w18;
1687                                const double tmp1_1 = A_01*w19;
1688                                const double tmp2_1 = A_10*w20;
1689                                const double tmp3_1 = A_11*w21;
1690                                const double tmp4_1 = A_00*w22;
1691                                const double tmp5_1 = tmp0_0*w19;
1692                                const double tmp6_1 = A_11*w23;
1693                                const double tmp7_1 = A_11*w25;
1694                                const double tmp8_1 = A_00*w24;
1695                                const double tmp9_1 = tmp0_0*w20;
1696                                const double tmp10_1 = A_01*w20;
1697                                const double tmp11_1 = A_11*w27;
1698                                const double tmp12_1 = A_00*w26;
1699                                const double tmp13_1 = A_10*w19;
1700                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1701                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1702                                EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1703                                EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1704                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1705                                EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1706                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1707                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1708                                EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1709                                EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1710                                EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1711                                EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1712                                EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1713                                EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1714                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1715                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1716                            }
1717                        }
1718                        ///////////////
1719                        // process B //
1720                        ///////////////
1721                        if (!B.isEmpty()) {
1722                            add_EM_S=true;
1723                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1724                            if (B.actsExpanded()) {
1725                                const double B_0_0 = B_p[INDEX2(0,0,2)];
1726                                const double B_1_0 = B_p[INDEX2(1,0,2)];
1727                                const double B_0_1 = B_p[INDEX2(0,1,2)];
1728                                const double B_1_1 = B_p[INDEX2(1,1,2)];
1729                                const double B_0_2 = B_p[INDEX2(0,2,2)];
1730                                const double B_1_2 = B_p[INDEX2(1,2,2)];
1731                                const double B_0_3 = B_p[INDEX2(0,3,2)];
1732                                const double B_1_3 = B_p[INDEX2(1,3,2)];
1733                                const double tmp0_0 = B_1_0 + B_1_1;
1734                                const double tmp1_0 = B_1_2 + B_1_3;
1735                                const double tmp2_0 = B_0_1 + B_0_3;
1736                                const double tmp3_0 = B_0_0 + B_0_2;
1737                                const double tmp63_1 = B_1_1*w42;
1738                                const double tmp79_1 = B_1_1*w40;
1739                                const double tmp37_1 = tmp3_0*w35;
1740                                const double tmp8_1 = tmp0_0*w32;
1741                                const double tmp71_1 = B_0_1*w34;
1742                                const double tmp19_1 = B_0_3*w31;
1743                                const double tmp15_1 = B_0_3*w34;
1744                                const double tmp9_1 = tmp3_0*w34;
1745                                const double tmp35_1 = B_1_0*w36;
1746                                const double tmp66_1 = B_0_3*w28;
1747                                const double tmp28_1 = B_1_0*w42;
1748                                const double tmp22_1 = B_1_0*w40;
1749                                const double tmp16_1 = B_1_2*w29;
1750                                const double tmp6_1 = tmp2_0*w35;
1751                                const double tmp55_1 = B_1_3*w40;
1752                                const double tmp50_1 = B_1_3*w42;
1753                                const double tmp7_1 = tmp1_0*w29;
1754                                const double tmp1_1 = tmp1_0*w32;
1755                                const double tmp57_1 = B_0_3*w30;
1756                                const double tmp18_1 = B_1_1*w32;
1757                                const double tmp53_1 = B_1_0*w41;
1758                                const double tmp61_1 = B_1_3*w36;
1759                                const double tmp27_1 = B_0_3*w38;
1760                                const double tmp64_1 = B_0_2*w30;
1761                                const double tmp76_1 = B_0_1*w38;
1762                                const double tmp39_1 = tmp2_0*w34;
1763                                const double tmp62_1 = B_0_1*w31;
1764                                const double tmp56_1 = B_0_0*w31;
1765                                const double tmp49_1 = B_1_1*w36;
1766                                const double tmp2_1 = B_0_2*w31;
1767                                const double tmp23_1 = B_0_2*w33;
1768                                const double tmp38_1 = B_1_1*w43;
1769                                const double tmp74_1 = B_1_2*w41;
1770                                const double tmp43_1 = B_1_1*w41;
1771                                const double tmp58_1 = B_0_2*w28;
1772                                const double tmp67_1 = B_0_0*w33;
1773                                const double tmp33_1 = tmp0_0*w39;
1774                                const double tmp4_1 = B_0_0*w28;
1775                                const double tmp20_1 = B_0_0*w30;
1776                                const double tmp13_1 = B_0_2*w38;
1777                                const double tmp65_1 = B_1_2*w43;
1778                                const double tmp0_1 = tmp0_0*w29;
1779                                const double tmp41_1 = tmp3_0*w33;
1780                                const double tmp73_1 = B_0_2*w37;
1781                                const double tmp69_1 = B_0_0*w38;
1782                                const double tmp48_1 = B_1_2*w39;
1783                                const double tmp59_1 = B_0_1*w33;
1784                                const double tmp17_1 = B_1_3*w41;
1785                                const double tmp5_1 = B_0_3*w33;
1786                                const double tmp3_1 = B_0_1*w30;
1787                                const double tmp21_1 = B_0_1*w28;
1788                                const double tmp42_1 = B_1_0*w29;
1789                                const double tmp54_1 = B_1_2*w32;
1790                                const double tmp60_1 = B_1_0*w39;
1791                                const double tmp32_1 = tmp1_0*w36;
1792                                const double tmp10_1 = B_0_1*w37;
1793                                const double tmp14_1 = B_0_0*w35;
1794                                const double tmp29_1 = B_0_1*w35;
1795                                const double tmp26_1 = B_1_2*w36;
1796                                const double tmp30_1 = B_1_3*w43;
1797                                const double tmp70_1 = B_0_2*w35;
1798                                const double tmp34_1 = B_1_3*w39;
1799                                const double tmp51_1 = B_1_0*w43;
1800                                const double tmp31_1 = B_0_2*w34;
1801                                const double tmp45_1 = tmp3_0*w28;
1802                                const double tmp11_1 = tmp1_0*w39;
1803                                const double tmp52_1 = B_1_1*w29;
1804                                const double tmp44_1 = B_1_3*w32;
1805                                const double tmp25_1 = B_1_1*w39;
1806                                const double tmp47_1 = tmp2_0*w33;
1807                                const double tmp72_1 = B_1_3*w29;
1808                                const double tmp40_1 = tmp2_0*w28;
1809                                const double tmp46_1 = B_1_2*w40;
1810                                const double tmp36_1 = B_1_2*w42;
1811                                const double tmp24_1 = B_0_0*w37;
1812                                const double tmp77_1 = B_0_3*w35;
1813                                const double tmp68_1 = B_0_3*w37;
1814                                const double tmp78_1 = B_0_0*w34;
1815                                const double tmp12_1 = tmp0_0*w36;
1816                                const double tmp75_1 = B_1_0*w32;
1817                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1818                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1819                                EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1820                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
1821                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1822                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1823                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1824                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1825                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1826                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1827                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1828                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1829                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
1830                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1831                                EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1832                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1833                            } else { // constant data
1834                                const double B_0 = B_p[0];
1835                                const double B_1 = B_p[1];
1836                                const double tmp0_1 = B_0*w44;
1837                                const double tmp1_1 = B_1*w45;
1838                                const double tmp2_1 = B_0*w46;
1839                                const double tmp3_1 = B_0*w47;
1840                                const double tmp4_1 = B_1*w48;
1841                                const double tmp5_1 = B_1*w49;
1842                                const double tmp6_1 = B_1*w50;
1843                                const double tmp7_1 = B_0*w51;
1844                                EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1845                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1846                                EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
1847                                EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
1848                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1849                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
1850                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
1851                                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
1852                                EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
1853                                EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
1854                                EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
1855                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1856                                EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
1857                                EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
1858                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
1859                                EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
1860                            }
1861                        }
1862                        ///////////////
1863                        // process C //
1864                        ///////////////
1865                        if (!C.isEmpty()) {
1866                            add_EM_S=true;
1867                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
1868                            if (C.actsExpanded()) {
1869                                const double C_0_0 = C_p[INDEX2(0,0,2)];
1870                                const double C_1_0 = C_p[INDEX2(1,0,2)];
1871                                const double C_0_1 = C_p[INDEX2(0,1,2)];
1872                                const double C_1_1 = C_p[INDEX2(1,1,2)];
1873                                const double C_0_2 = C_p[INDEX2(0,2,2)];
1874                                const double C_1_2 = C_p[INDEX2(1,2,2)];
1875                                const double C_0_3 = C_p[INDEX2(0,3,2)];
1876                                const double C_1_3 = C_p[INDEX2(1,3,2)];
1877                                const double tmp0_0 = C_1_0 + C_1_1;
1878                                const double tmp1_0 = C_1_2 + C_1_3;
1879                                const double tmp2_0 = C_0_1 + C_0_3;
1880                                const double tmp3_0 = C_0_0 + C_0_2;
1881                                const double tmp64_1 = C_0_2*w30;
1882                                const double tmp14_1 = C_0_2*w28;
1883                                const double tmp19_1 = C_0_3*w31;
1884                                const double tmp22_1 = C_1_0*w40;
1885                                const double tmp37_1 = tmp3_0*w35;
1886                                const double tmp29_1 = C_0_1*w35;
1887                                const double tmp73_1 = C_0_2*w37;
1888                                const double tmp74_1 = C_1_2*w41;
1889                                const double tmp52_1 = C_1_3*w39;
1890                                const double tmp25_1 = C_1_1*w39;
1891                                const double tmp62_1 = C_0_1*w31;
1892                                const double tmp79_1 = C_1_1*w40;
1893                                const double tmp43_1 = C_1_1*w36;
1894                                const double tmp27_1 = C_0_3*w38;
1895                                const double tmp28_1 = C_1_0*w42;
1896                                const double tmp63_1 = C_1_1*w42;
1897                                const double tmp59_1 = C_0_3*w34;
1898                                const double tmp72_1 = C_1_3*w29;
1899                                const double tmp40_1 = tmp2_0*w35;
1900                                const double tmp13_1 = C_0_3*w30;
1901                                const double tmp51_1 = C_1_2*w40;
1902                                const double tmp54_1 = C_1_2*w42;
1903                                const double tmp12_1 = C_0_0*w31;
1904                                const double tmp2_1 = tmp1_0*w32;
1905                                const double tmp68_1 = C_0_2*w31;
1906                                const double tmp75_1 = C_1_0*w32;
1907                                const double tmp49_1 = C_1_1*w41;
1908                                const double tmp4_1 = C_0_2*w35;
1909                                const double tmp66_1 = C_0_3*w28;
1910                                const double tmp56_1 = C_0_1*w37;
1911                                const double tmp5_1 = C_0_1*w34;
1912                                const double tmp38_1 = tmp2_0*w34;
1913                                const double tmp76_1 = C_0_1*w38;
1914                                const double tmp21_1 = C_0_1*w28;
1915                                const double tmp69_1 = C_0_1*w30;
1916                                const double tmp53_1 = C_1_0*w36;
1917                                const double tmp42_1 = C_1_2*w39;
1918                                const double tmp32_1 = tmp1_0*w29;
1919                                const double tmp45_1 = C_1_0*w43;
1920                                const double tmp33_1 = tmp0_0*w32;
1921                                const double tmp35_1 = C_1_0*w41;
1922                                const double tmp26_1 = C_1_2*w36;
1923                                const double tmp67_1 = C_0_0*w33;
1924                                const double tmp31_1 = C_0_2*w34;
1925                                const double tmp20_1 = C_0_0*w30;
1926                                const double tmp70_1 = C_0_0*w28;
1927                                const double tmp8_1 = tmp0_0*w39;
1928                                const double tmp30_1 = C_1_3*w43;
1929                                const double tmp0_1 = tmp0_0*w29;
1930                                const double tmp17_1 = C_1_3*w41;
1931                                const double tmp58_1 = C_0_0*w35;
1932                                const double tmp9_1 = tmp3_0*w33;
1933                                const double tmp61_1 = C_1_3*w36;
1934                                const double tmp41_1 = tmp3_0*w34;
1935                                const double tmp50_1 = C_1_3*w32;
1936                                const double tmp18_1 = C_1_1*w32;
1937                                const double tmp6_1 = tmp1_0*w36;
1938                                const double tmp3_1 = C_0_0*w38;
1939                                const double tmp34_1 = C_1_1*w29;
1940                                const double tmp77_1 = C_0_3*w35;
1941                                const double tmp65_1 = C_1_2*w43;
1942                                const double tmp71_1 = C_0_3*w33;
1943                                const double tmp55_1 = C_1_1*w43;
1944                                const double tmp46_1 = tmp3_0*w28;
1945                                const double tmp24_1 = C_0_0*w37;
1946                                const double tmp10_1 = tmp1_0*w39;
1947                                const double tmp48_1 = C_1_0*w29;
1948                                const double tmp15_1 = C_0_1*w33;
1949                                const double tmp36_1 = C_1_2*w32;
1950                                const double tmp60_1 = C_1_0*w39;
1951                                const double tmp47_1 = tmp2_0*w33;
1952                                const double tmp16_1 = C_1_2*w29;
1953                                const double tmp1_1 = C_0_3*w37;
1954                                const double tmp7_1 = tmp2_0*w28;
1955                                const double tmp39_1 = C_1_3*w40;
1956                                const double tmp44_1 = C_1_3*w42;
1957                                const double tmp57_1 = C_0_2*w38;
1958                                const double tmp78_1 = C_0_0*w34;
1959                                const double tmp11_1 = tmp0_0*w36;
1960                                const double tmp23_1 = C_0_2*w33;
1961                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1962                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1963                                EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1964                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
1965                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1966                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1967                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1968                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1969                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1970                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1971                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1972                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1973                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
1974                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1975                                EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1976                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1977                            } else { // constant data
1978                                const double C_0 = C_p[0];
1979                                const double C_1 = C_p[1];
1980                                const double tmp0_1 = C_0*w47;
1981                                const double tmp1_1 = C_1*w45;
1982                                const double tmp2_1 = C_1*w48;
1983                                const double tmp3_1 = C_0*w51;
1984                                const double tmp4_1 = C_0*w44;
1985                                const double tmp5_1 = C_1*w49;
1986                                const double tmp6_1 = C_1*w50;
1987                                const double tmp7_1 = C_0*w46;
1988                                EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1989                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
1990                                EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
1991                                EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
1992                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1993                                EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
1994                                EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
1995                                EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
1996                                EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
1997                                EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
1998                                EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
1999                                EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
2000                                EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
2001                                EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
2002                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2003                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
2004                            }
2005                        }
2006                        ///////////////
2007                        // process D //
2008                        ///////////////
2009                        if (!D.isEmpty()) {
2010                            add_EM_S=true;
2011                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2012                            if (D.actsExpanded()) {
2013                                const double D_0 = D_p[0];
2014                                const double D_1 = D_p[1];
2015                                const double D_2 = D_p[2];
2016                                const double D_3 = D_p[3];
2017                                const double tmp0_0 = D_0 + D_1;
2018                                const double tmp1_0 = D_2 + D_3;
2019                                const double tmp2_0 = D_0 + D_1 + D_2 + D_3;
2020                                const double tmp3_0 = D_1 + D_2;
2021                                const double tmp4_0 = D_1 + D_3;
2022                                const double tmp5_0 = D_0 + D_2;
2023                                const double tmp6_0 = D_0 + D_3;
2024                                const double tmp0_1 = tmp0_0*w52;
2025                                const double tmp1_1 = tmp1_0*w53;
2026                                const double tmp2_1 = tmp2_0*w54;
2027                                const double tmp3_1 = tmp1_0*w52;
2028                                const double tmp4_1 = tmp0_0*w53;
2029                                const double tmp5_1 = tmp3_0*w54;
2030                                const double tmp6_1 = D_0*w55;
2031                                const double tmp7_1 = D_3*w56;
2032                                const double tmp8_1 = D_3*w55;
2033                                const double tmp9_1 = D_0*w56;
2034                                const double tmp10_1 = tmp4_0*w52;
2035                                const double tmp11_1 = tmp5_0*w53;
2036                                const double tmp12_1 = tmp5_0*w52;
2037                                const double tmp13_1 = tmp4_0*w53;
2038                                const double tmp14_1 = tmp6_0*w54;
2039                                const double tmp15_1 = D_2*w55;
2040                                const double tmp16_1 = D_1*w56;
2041                                const double tmp17_1 = D_1*w55;
2042                                const double tmp18_1 = D_2*w56;
2043                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2044                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2045                                EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
2046                                EM_S[INDEX2(3,0,4)]+=tmp2_1;
2047                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2048                                EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
2049                                EM_S[INDEX2(2,1,4)]+=tmp2_1;
2050                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
2051                                EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
2052                                EM_S[INDEX2(1,2,4)]+=tmp2_1;
2053                                EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
2054                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
2055                                EM_S[INDEX2(0,3,4)]+=tmp2_1;
2056                                EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
2057                                EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2058                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2059                            } else { // constant data
2060                                const double tmp0_1 = D_p[0]*w57;
2061                                const double tmp1_1 = D_p[0]*w58;
2062                                const double tmp2_1 = D_p[0]*w59;
2063                                EM_S[INDEX2(0,0,4)]+=tmp2_1;
2064                                EM_S[INDEX2(1,0,4)]+=tmp0_1;
2065                                EM_S[INDEX2(2,0,4)]+=tmp0_1;
2066                                EM_S[INDEX2(3,0,4)]+=tmp1_1;
2067                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
2068                                EM_S[INDEX2(1,1,4)]+=tmp2_1;
2069                                EM_S[INDEX2(2,1,4)]+=tmp1_1;
2070                                EM_S[INDEX2(3,1,4)]+=tmp0_1;
2071                                EM_S[INDEX2(0,2,4)]+=tmp0_1;
2072                                EM_S[INDEX2(1,2,4)]+=tmp1_1;
2073                                EM_S[INDEX2(2,2,4)]+=tmp2_1;
2074                                EM_S[INDEX2(3,2,4)]+=tmp0_1;
2075                                EM_S[INDEX2(0,3,4)]+=tmp1_1;
2076                                EM_S[INDEX2(1,3,4)]+=tmp0_1;
2077                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
2078                                EM_S[INDEX2(3,3,4)]+=tmp2_1;
2079                            }
2080                        }
2081                        ///////////////
2082                        // process X //
2083                        ///////////////
2084                        if (!X.isEmpty()) {
2085                            add_EM_F=true;
2086                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2087                            if (X.actsExpanded()) {
2088                                const double X_0_0 = X_p[INDEX2(0,0,2)];
2089                                const double X_1_0 = X_p[INDEX2(1,0,2)];
2090                                const double X_0_1 = X_p[INDEX2(0,1,2)];
2091                                const double X_1_1 = X_p[INDEX2(1,1,2)];
2092                                const double X_0_2 = X_p[INDEX2(0,2,2)];
2093                                const double X_1_2 = X_p[INDEX2(1,2,2)];
2094                                const double X_0_3 = X_p[INDEX2(0,3,2)];
2095                                const double X_1_3 = X_p[INDEX2(1,3,2)];
2096                                const double tmp0_0 = X_0_2 + X_0_3;
2097                                const double tmp1_0 = X_1_1 + X_1_3;
2098                                const double tmp2_0 = X_1_0 + X_1_2;
2099                                const double tmp3_0 = X_0_0 + X_0_1;
2100                                const double tmp0_1 = tmp0_0*w63;
2101                                const double tmp1_1 = tmp1_0*w62;
2102                                const double tmp2_1 = tmp2_0*w61;
2103                                const double tmp3_1 = tmp3_0*w60;
2104                                const double tmp4_1 = tmp0_0*w65;
2105                                const double tmp5_1 = tmp3_0*w64;
2106                                const double tmp6_1 = tmp2_0*w62;
2107                                const double tmp7_1 = tmp1_0*w61;
2108                                const double tmp8_1 = tmp2_0*w66;
2109                                const double tmp9_1 = tmp3_0*w63;
2110                                const double tmp10_1 = tmp0_0*w60;
2111                                const double tmp11_1 = tmp1_0*w67;
2112                                const double tmp12_1 = tmp1_0*w66;
2113                                const double tmp13_1 = tmp3_0*w65;
2114                                const double tmp14_1 = tmp0_0*w64;
2115                                const double tmp15_1 = tmp2_0*w67;
2116                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2117                                EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2118                                EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2119                                EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2120                            } else { // constant data
2121                                const double tmp0_1 = X_p[1]*w69;
2122                                const double tmp1_1 = X_p[0]*w68;
2123                                const double tmp2_1 = X_p[0]*w70;
2124                                const double tmp3_1 = X_p[1]*w71;
2125                                EM_F[0]+=tmp0_1 + tmp1_1;
2126                                EM_F[1]+=tmp0_1 + tmp2_1;
2127                                EM_F[2]+=tmp1_1 + tmp3_1;
2128                                EM_F[3]+=tmp2_1 + tmp3_1;
2129                            }
2130                        }
2131                        ///////////////
2132                        // process Y //
2133                        ///////////////
2134                        if (!Y.isEmpty()) {
2135                            add_EM_F=true;
2136                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2137                            if (Y.actsExpanded()) {
2138                                const double Y_0 = Y_p[0];
2139                                const double Y_1 = Y_p[1];
2140                                const double Y_2 = Y_p[2];
2141                                const double Y_3 = Y_p[3];
2142                                const double tmp0_0 = Y_1 + Y_2;
2143                                const double tmp1_0 = Y_0 + Y_3;
2144                                const double tmp0_1 = Y_0*w72;
2145                                const double tmp1_1 = tmp0_0*w73;
2146                                const double tmp2_1 = Y_3*w74;
2147                                const double tmp3_1 = Y_1*w72;
2148                                const double tmp4_1 = tmp1_0*w73;
2149                                const double tmp5_1 = Y_2*w74;
2150                                const double tmp6_1 = Y_2*w72;
2151                                const double tmp7_1 = Y_1*w74;
2152                                const double tmp8_1 = Y_3*w72;
2153                                const double tmp9_1 = Y_0*w74;
2154                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2155                                EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2156                                EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2157                                EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2158                            } else { // constant data
2159                                const double tmp0_1 = Y_p[0]*w75;
2160                                EM_F[0]+=tmp0_1;
2161                                EM_F[1]+=tmp0_1;
2162                                EM_F[2]+=tmp0_1;
2163                                EM_F[3]+=tmp0_1;
2164                            }
2165                        }
2166    
2167                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2168                        const index_t firstNode=m_N0*k1+k0;
2169                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2170                    } // end k0 loop
2171                } // end k1 loop
2172            } // end of colouring
2173        } // end of parallel region
2174    }
2175    
2176    //protected
2177    void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2178            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2179            const escript::Data& C, const escript::Data& D,
2180            const escript::Data& X, const escript::Data& Y) const
2181    {
2182        const double h0 = m_l0/m_gNE0;
2183        const double h1 = m_l1/m_gNE1;
2184        const double w0 = -.25*h1/h0;
2185        const double w1 = .25;
2186        const double w2 = -.25;
2187        const double w3 = .25*h0/h1;
2188        const double w4 = -.25*h0/h1;
2189        const double w5 = .25*h1/h0;
2190        const double w6 = -.125*h1;
2191        const double w7 = -.125*h0;
2192        const double w8 = .125*h1;
2193        const double w9 = .125*h0;
2194        const double w10 = .0625*h0*h1;
2195        const double w11 = -.5*h1;
2196        const double w12 = -.5*h0;
2197        const double w13 = .5*h1;
2198        const double w14 = .5*h0;
2199        const double w15 = .25*h0*h1;
2200    
2201        rhs.requireWrite();
2202    #pragma omp parallel
2203        {
2204            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2205    #pragma omp for
2206                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2207                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2208                        bool add_EM_S=false;
2209                        bool add_EM_F=false;
2210                        vector<double> EM_S(4*4, 0);
2211                        vector<double> EM_F(4, 0);
2212                        const index_t e = k0 + m_NE0*k1;
2213                        ///////////////
2214                        // process A //
2215                        ///////////////
2216                        if (!A.isEmpty()) {
2217                            add_EM_S=true;
2218                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2219                            const double A_00 = A_p[INDEX2(0,0,2)];
2220                            const double A_10 = A_p[INDEX2(1,0,2)];
2221                            const double A_01 = A_p[INDEX2(0,1,2)];
2222                            const double A_11 = A_p[INDEX2(1,1,2)];
2223                            const double tmp0_0 = A_01 + A_10;
2224                            const double tmp0_1 = A_11*w3;
2225                            const double tmp1_1 = A_00*w0;
2226                            const double tmp2_1 = A_01*w1;
2227                            const double tmp3_1 = A_10*w2;
2228                            const double tmp4_1 = tmp0_0*w1;
2229                            const double tmp5_1 = A_11*w4;
2230                            const double tmp6_1 = A_00*w5;
2231                            const double tmp7_1 = tmp0_0*w2;
2232                            const double tmp8_1 = A_10*w1;
2233                            const double tmp9_1 = A_01*w2;
2234                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2235                            EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2236                            EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2237                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2238                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2239                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2240                            EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2241                            EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2242                            EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2243                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2244                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2245                            EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2246                            EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2247                            EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2248                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2249                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2250                        }
2251                        ///////////////
2252                        // process B //
2253                        ///////////////
2254                        if (!B.isEmpty()) {
2255                            add_EM_S=true;
2256                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2257                            const double tmp2_1 = B_p[0]*w8;
2258                            const double tmp0_1 = B_p[0]*w6;
2259                            const double tmp3_1 = B_p[1]*w9;
2260                            const double tmp1_1 = B_p[1]*w7;
2261                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2262                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2263                            EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
2264                            EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;
2265                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2266                            EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;
2267                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;
2268                            EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;
2269                            EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;
2270                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
2271                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;
2272                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2273                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;
2274                            EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;
2275                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;
2276                            EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
2277                        }
2278                        ///////////////
2279                        // process C //
2280                        ///////////////
2281                        if (!C.isEmpty()) {
2282                            add_EM_S=true;
2283                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2284                            const double tmp1_1 = C_p[1]*w7;
2285                            const double tmp0_1 = C_p[0]*w8;
2286                            const double tmp3_1 = C_p[0]*w6;
2287                            const double tmp2_1 = C_p[1]*w9;
2288                            EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2289                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2290                            EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
2291                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
2292                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2293                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;
2294                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;
2295                            EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;
2296                            EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;
2297                            EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
2298                            EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
2299                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2300                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;
2301                            EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;
2302                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2303                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;
2304                        }
2305                        ///////////////
2306                        // process D //
2307                        ///////////////
2308                        if (!D.isEmpty()) {
2309                            add_EM_S=true;
2310                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2311                            const double tmp0_1 = D_p[0]*w10;
2312                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
2313                            EM_S[INDEX2(1,0,4)]+=tmp0_1;
2314                            EM_S[INDEX2(2,0,4)]+=tmp0_1;
2315                            EM_S[INDEX2(3,0,4)]+=tmp0_1;
2316                            EM_S[INDEX2(0,1,4)]+=tmp0_1;
2317                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
2318                            EM_S[INDEX2(2,1,4)]+=tmp0_1;
2319                            EM_S[INDEX2(3,1,4)]+=tmp0_1;
2320                            EM_S[INDEX2(0,2,4)]+=tmp0_1;
2321                            EM_S[INDEX2(1,2,4)]+=tmp0_1;
2322                            EM_S[INDEX2(2,2,4)]+=tmp0_1;
2323                            EM_S[INDEX2(3,2,4)]+=tmp0_1;
2324                            EM_S[INDEX2(0,3,4)]+=tmp0_1;
2325                            EM_S[INDEX2(1,3,4)]+=tmp0_1;
2326                            EM_S[INDEX2(2,3,4)]+=tmp0_1;
2327                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
2328                        }
2329                        ///////////////
2330                        // process X //
2331                        ///////////////
2332                        if (!X.isEmpty()) {
2333                            add_EM_F=true;
2334                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2335                            const double tmp0_1 = X_p[0]*w11;
2336                            const double tmp2_1 = X_p[0]*w13;
2337                            const double tmp1_1 = X_p[1]*w12;
2338                            const double tmp3_1 = X_p[1]*w14;
2339                            EM_F[0]+=tmp0_1 + tmp1_1;
2340                            EM_F[1]+=tmp1_1 + tmp2_1;
2341                            EM_F[2]+=tmp0_1 + tmp3_1;
2342                            EM_F[3]+=tmp2_1 + tmp3_1;
2343                        }
2344                        ///////////////
2345                        // process Y //
2346                        ///////////////
2347                        if (!Y.isEmpty()) {
2348                            add_EM_F=true;
2349                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2350                            const double tmp0_1 = Y_p[0]*w15;
2351                            EM_F[0]+=tmp0_1;
2352                            EM_F[1]+=tmp0_1;
2353                            EM_F[2]+=tmp0_1;
2354                            EM_F[3]+=tmp0_1;
2355                        }
2356    
2357                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2358                        const index_t firstNode=m_N0*k1+k0;
2359                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2360                    } // end k0 loop
2361                } // end k1 loop
2362            } // end of colouring
2363        } // end of parallel region
2364    }
2365    
2366    //protected
2367    void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2368            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2369            const escript::Data& C, const escript::Data& D,
2370            const escript::Data& X, const escript::Data& Y) const
2371    {
2372        const double h0 = m_l0/m_gNE0;
2373        const double h1 = m_l1/m_gNE1;
2374        dim_t numEq, numComp;
2375        if (!mat)
2376            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
2377        else {
2378            numEq=mat->logical_row_block_size;
2379            numComp=mat->logical_col_block_size;
2380        }
2381    
2382        const double w0 = -0.1555021169820365539*h1/h0;
2383        const double w1 = 0.041666666666666666667;
2384        const double w2 = -0.15550211698203655390;
2385        const double w3 = 0.041666666666666666667*h0/h1;
2386        const double w4 = 0.15550211698203655390;
2387        const double w5 = -0.041666666666666666667;
2388        const double w6 = -0.01116454968463011277*h1/h0;
2389        const double w7 = 0.011164549684630112770;
2390        const double w8 = -0.011164549684630112770;
2391        const double w9 = -0.041666666666666666667*h1/h0;
2392        const double w10 = -0.041666666666666666667*h0/h1;
2393        const double w11 = 0.1555021169820365539*h1/h0;
2394        const double w12 = 0.1555021169820365539*h0/h1;
2395        const double w13 = 0.01116454968463011277*h0/h1;
2396        const double w14 = 0.01116454968463011277*h1/h0;
2397        const double w15 = 0.041666666666666666667*h1/h0;
2398        const double w16 = -0.01116454968463011277*h0/h1;
2399        const double w17 = -0.1555021169820365539*h0/h1;
2400        const double w18 = -0.33333333333333333333*h1/h0;
2401        const double w19 = 0.25000000000000000000;
2402        const double w20 = -0.25000000000000000000;
2403        const double w21 = 0.16666666666666666667*h0/h1;
2404        const double w22 = -0.16666666666666666667*h1/h0;
2405        const double w23 = -0.16666666666666666667*h0/h1;
2406        const double w24 = 0.33333333333333333333*h1/h0;
2407        const double w25 = 0.33333333333333333333*h0/h1;
2408        const double w26 = 0.16666666666666666667*h1/h0;
2409        const double w27 = -0.33333333333333333333*h0/h1;
2410        const double w28 = -0.032861463941450536761*h1;
2411        const double w29 = -0.032861463941450536761*h0;
2412        const double w30 = -0.12264065304058601714*h1;
2413        const double w31 = -0.0023593469594139828636*h1;
2414        const double w32 = -0.008805202725216129906*h0;
2415        const double w33 = -0.008805202725216129906*h1;
2416        const double w34 = 0.032861463941450536761*h1;
2417        const double w35 = 0.008805202725216129906*h1;
2418        const double w36 = 0.008805202725216129906*h0;
2419        const double w37 = 0.0023593469594139828636*h1;
2420        const double w38 = 0.12264065304058601714*h1;
2421        const double w39 = 0.032861463941450536761*h0;
2422        const double w40 = -0.12264065304058601714*h0;
2423        const double w41 = -0.0023593469594139828636*h0;
2424        const double w42 = 0.0023593469594139828636*h0;
2425        const double w43 = 0.12264065304058601714*h0;
2426        const double w44 = -0.16666666666666666667*h1;
2427        const double w45 = -0.083333333333333333333*h0;
2428        const double w46 = 0.083333333333333333333*h1;
2429        const double w47 = 0.16666666666666666667*h1;
2430        const double w48 = 0.083333333333333333333*h0;
2431        const double w49 = -0.16666666666666666667*h0;
2432        const double w50 = 0.16666666666666666667*h0;
2433        const double w51 = -0.083333333333333333333*h1;
2434        const double w52 = 0.025917019497006092316*h0*h1;
2435        const double w53 = 0.0018607582807716854616*h0*h1;
2436        const double w54 = 0.0069444444444444444444*h0*h1;
2437        const double w55 = 0.09672363354357992482*h0*h1;
2438        const double w56 = 0.00049858867864229740201*h0*h1;
2439        const double w57 = 0.055555555555555555556*h0*h1;
2440        const double w58 = 0.027777777777777777778*h0*h1;
2441        const double w59 = 0.11111111111111111111*h0*h1;
2442        const double w60 = -0.19716878364870322056*h1;
2443        const double w61 = -0.19716878364870322056*h0;
2444        const double w62 = -0.052831216351296779436*h0;
2445        const double w63 = -0.052831216351296779436*h1;
2446        const double w64 = 0.19716878364870322056*h1;
2447        const double w65 = 0.052831216351296779436*h1;
2448        const double w66 = 0.19716878364870322056*h0;
2449        const double w67 = 0.052831216351296779436*h0;
2450        const double w68 = -0.5*h1;
2451        const double w69 = -0.5*h0;
2452        const double w70 = 0.5*h1;
2453        const double w71 = 0.5*h0;
2454        const double w72 = 0.1555021169820365539*h0*h1;
2455        const double w73 = 0.041666666666666666667*h0*h1;
2456        const double w74 = 0.01116454968463011277*h0*h1;
2457        const double w75 = 0.25*h0*h1;
2458    
2459        rhs.requireWrite();
2460    #pragma omp parallel
2461        {
2462            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2463    #pragma omp for
2464                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2465                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2466                        bool add_EM_S=false;
2467                        bool add_EM_F=false;
2468                        vector<double> EM_S(4*4*numEq*numComp, 0);
2469                        vector<double> EM_F(4*numEq, 0);
2470                        const index_t e = k0 + m_NE0*k1;
2471                        ///////////////
2472                        // process A //
2473                        ///////////////
2474                        if (!A.isEmpty()) {
2475                            add_EM_S=true;
2476                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2477                            if (A.actsExpanded()) {
2478                                for (index_t k=0; k<numEq; k++) {
2479                                    for (index_t m=0; m<numComp; m++) {
2480                                        const double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];
2481                                        const double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];
2482                                        const double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];
2483                                        const double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];
2484                                        const double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];
2485                                        const double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];
2486                                        const double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];
2487                                        const double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];
2488                                        const double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];
2489                                        const double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];
2490                                        const double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];
2491                                        const double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];
2492                                        const double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];
2493                                        const double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];
2494                                        const double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];
2495                                        const double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];
2496                                        const double tmp0_0 = A_01_0 + A_01_3;
2497                                        const double tmp1_0 = A_00_0 + A_00_1;
2498                                        const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
2499                                        const double tmp3_0 = A_00_2 + A_00_3;
2500                                        const double tmp4_0 = A_10_1 + A_10_2;
2501                                        const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
2502                                        const double tmp6_0 = A_01_3 + A_10_0;
2503                                        const double tmp7_0 = A_01_0 + A_10_3;
2504                                        const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
2505                                        const double tmp9_0 = A_01_0 + A_10_0;
2506                                        const double tmp10_0 = A_01_3 + A_10_3;
2507                                        const double tmp11_0 = A_11_1 + A_11_3;
2508                                        const double tmp12_0 = A_11_0 + A_11_2;
2509                                        const double tmp13_0 = A_01_2 + A_10_1;
2510                                        const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
2511                                        const double tmp15_0 = A_01_1 + A_10_2;
2512                                        const double tmp16_0 = A_10_0 + A_10_3;
2513                                        const double tmp17_0 = A_01_1 + A_01_2;
2514                                        const double tmp18_0 = A_01_1 + A_10_1;
2515                                        const double tmp19_0 = A_01_2 + A_10_2;
2516                                        const double tmp0_1 = A_10_3*w8;
2517                                        const double tmp1_1 = tmp0_0*w1;
2518                                        const double tmp2_1 = A_01_1*w4;
2519                                        const double tmp3_1 = tmp1_0*w0;
2520                                        const double tmp4_1 = A_01_2*w7;
2521                                        const double tmp5_1 = tmp2_0*w3;
2522                                        const double tmp6_1 = tmp3_0*w6;
2523                                        const double tmp7_1 = A_10_0*w2;
2524                                        const double tmp8_1 = tmp4_0*w5;
2525                                        const double tmp9_1 = tmp2_0*w10;
2526                                        const double tmp10_1 = tmp5_0*w9;
2527                                        const double tmp11_1 = tmp6_0*w7;
2528                                        const double tmp12_1 = tmp7_0*w4;
2529                                        const double tmp13_1 = tmp8_0*w1;
2530                                        const double tmp14_1 = A_10_0*w8;
2531                                        const double tmp15_1 = A_01_2*w4;
2532                                        const double tmp16_1 = tmp3_0*w0;
2533                                        const double tmp17_1 = A_01_1*w7;
2534                                        const double tmp18_1 = tmp1_0*w6;
2535                                        const double tmp19_1 = A_10_3*w2;
2536                                        const double tmp20_1 = tmp9_0*w4;
2537                                        const double tmp21_1 = tmp1_0*w11;
2538                                        const double tmp22_1 = tmp10_0*w7;
2539                                        const double tmp23_1 = tmp3_0*w14;
2540                                        const double tmp24_1 = tmp11_0*w13;
2541                                        const double tmp25_1 = tmp12_0*w12;
2542                                        const double tmp26_1 = tmp10_0*w4;
2543                                        const double tmp27_1 = tmp3_0*w11;
2544                                        const double tmp28_1 = tmp9_0*w7;
2545                                        const double tmp29_1 = tmp1_0*w14;
2546                                        const double tmp30_1 = tmp12_0*w13;
2547                                        const double tmp31_1 = tmp11_0*w12;
2548                                        const double tmp32_1 = tmp13_0*w2;
2549                                        const double tmp33_1 = tmp14_0*w5;
2550                                        const double tmp34_1 = tmp15_0*w8;
2551                                        const double tmp35_1 = A_01_0*w8;
2552                                        const double tmp36_1 = tmp16_0*w1;
2553                                        const double tmp37_1 = A_10_1*w4;
2554                                        const double tmp38_1 = tmp5_0*w15;
2555                                        const double tmp39_1 = A_10_2*w7;
2556                                        const double tmp40_1 = tmp11_0*w17;
2557                                        const double tmp41_1 = A_01_3*w2;
2558                                        const double tmp42_1 = tmp12_0*w16;
2559                                        const double tmp43_1 = tmp17_0*w5;
2560                                        const double tmp44_1 = tmp7_0*w7;
2561                                        const double tmp45_1 = tmp6_0*w4;
2562                                        const double tmp46_1 = A_01_3*w8;
2563                                        const double tmp47_1 = A_10_2*w4;
2564                                        const double tmp48_1 = A_10_1*w7;
2565                                        const double tmp49_1 = tmp12_0*w17;
2566                                        const double tmp50_1 = A_01_0*w2;
2567                                        const double tmp51_1 = tmp11_0*w16;
2568                                        const double tmp52_1 = tmp18_0*w8;
2569                                        const double tmp53_1 = tmp19_0*w2;
2570                                        const double tmp54_1 = tmp13_0*w8;
2571                                        const double tmp55_1 = tmp15_0*w2;
2572                                        const double tmp56_1 = tmp19_0*w8;
2573                                        const double tmp57_1 = tmp18_0*w2;
2574                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
2575                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
2576                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
2577                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
2578                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
2579                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
2580                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
2581                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
2582                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2583                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
2584                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
2585                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
2586                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
2587                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
2588                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
2589                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2590                                    }
2591                                }
2592                            } else { // constant data
2593                                for (index_t k=0; k<numEq; k++) {
2594                                    for (index_t m=0; m<numComp; m++) {
2595                                        const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
2596                                        const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
2597                                        const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
2598                                        const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
2599                                        const double tmp0_0 = A_01 + A_10;
2600                                        const double tmp0_1 = A_00*w18;
2601                                        const double tmp1_1 = A_01*w19;
2602                                        const double tmp2_1 = A_10*w20;
2603                                        const double tmp3_1 = A_11*w21;
2604                                        const double tmp4_1 = A_00*w22;
2605                                        const double tmp5_1 = tmp0_0*w19;
2606                                        const double tmp6_1 = A_11*w23;
2607                                        const double tmp7_1 = A_11*w25;
2608                                        const double tmp8_1 = A_00*w24;
2609                                        const double tmp9_1 = tmp0_0*w20;
2610                                        const double tmp10_1 = A_01*w20;
2611                                        const double tmp11_1 = A_11*w27;
2612                                        const double tmp12_1 = A_00*w26;
2613                                        const double tmp13_1 = A_10*w19;
2614                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2615                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2616                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
2617                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2618                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2619                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2620                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
2621                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
2622                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
2623                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
2624                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2625                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2626                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2627                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
2628                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2629                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2630                                    }
2631                                }
2632                            }
2633                        }
2634                        ///////////////
2635                        // process B //
2636                        ///////////////
2637                        if (!B.isEmpty()) {
2638                            add_EM_S=true;
2639                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2640                            if (B.actsExpanded()) {
2641                                for (index_t k=0; k<numEq; k++) {
2642                                    for (index_t m=0; m<numComp; m++) {
2643                                        const double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
2644                                        const double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
2645                                        const double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
2646                                        const double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
2647                                        const double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];
2648                                        const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
2649                                        const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
2650                                        const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
2651                                        const double tmp0_0 = B_1_0 + B_1_1;
2652                                        const double tmp1_0 = B_1_2 + B_1_3;
2653                                        const double tmp2_0 = B_0_1 + B_0_3;
2654                                        const double tmp3_0 = B_0_0 + B_0_2;
2655                                        const double tmp63_1 = B_1_1*w42;
2656                                        const double tmp79_1 = B_1_1*w40;
2657                                        const double tmp37_1 = tmp3_0*w35;
2658                                        const double tmp8_1 = tmp0_0*w32;
2659                                        const double tmp71_1 = B_0_1*w34;
2660                                        const double tmp19_1 = B_0_3*w31;
2661                                        const double tmp15_1 = B_0_3*w34;
2662                                        const double tmp9_1 = tmp3_0*w34;
2663                                        const double tmp35_1 = B_1_0*w36;
2664                                        const double tmp66_1 = B_0_3*w28;
2665                                        const double tmp28_1 = B_1_0*w42;
2666                                        const double tmp22_1 = B_1_0*w40;
2667                                        const double tmp16_1 = B_1_2*w29;
2668                                        const double tmp6_1 = tmp2_0*w35;
2669                                        const double tmp55_1 = B_1_3*w40;
2670                                        const double tmp50_1 = B_1_3*w42;
2671                                        const double tmp7_1 = tmp1_0*w29;
2672                                        const double tmp1_1 = tmp1_0*w32;
2673                                        const double tmp57_1 = B_0_3*w30;
2674                                        const double tmp18_1 = B_1_1*w32;
2675                                        const double tmp53_1 = B_1_0*w41;
2676                                        const double tmp61_1 = B_1_3*w36;
2677                                        const double tmp27_1 = B_0_3*w38;
2678                                        const double tmp64_1 = B_0_2*w30;
2679                                        const double tmp76_1 = B_0_1*w38;
2680                                        const double tmp39_1 = tmp2_0*w34;
2681                                        const double tmp62_1 = B_0_1*w31;
2682                                        const double tmp56_1 = B_0_0*w31;
2683                                        const double tmp49_1 = B_1_1*w36;
2684                                        const double tmp2_1 = B_0_2*w31;
2685                                        const double tmp23_1 = B_0_2*w33;
2686                                        const double tmp38_1 = B_1_1*w43;
2687                                        const double tmp74_1 = B_1_2*w41;
2688                                        const double tmp43_1 = B_1_1*w41;
2689                                        const double tmp58_1 = B_0_2*w28;
2690                                        const double tmp67_1 = B_0_0*w33;
2691                                        const double tmp33_1 = tmp0_0*w39;
2692                                        const double tmp4_1 = B_0_0*w28;
2693                                        const double tmp20_1 = B_0_0*w30;
2694                                        const double tmp13_1 = B_0_2*w38;
2695                                        const double tmp65_1 = B_1_2*w43;
2696                                        const double tmp0_1 = tmp0_0*w29;
2697                                        const double tmp41_1 = tmp3_0*w33;
2698                                        const double tmp73_1 = B_0_2*w37;
2699                                        const double tmp69_1 = B_0_0*w38;
2700                                        const double tmp48_1 = B_1_2*w39;
2701                                        const double tmp59_1 = B_0_1*w33;
2702                                        const double tmp17_1 = B_1_3*w41;
2703                                        const double tmp5_1 = B_0_3*w33;
2704                                        const double tmp3_1 = B_0_1*w30;
2705                                        const double tmp21_1 = B_0_1*w28;
2706                                        const double tmp42_1 = B_1_0*w29;
2707                                        const double tmp54_1 = B_1_2*w32;
2708                                        const double tmp60_1 = B_1_0*w39;
2709                                        const double tmp32_1 = tmp1_0*w36;
2710                                        const double tmp10_1 = B_0_1*w37;
2711                                        const double tmp14_1 = B_0_0*w35;
2712                                        const double tmp29_1 = B_0_1*w35;
2713                                        const double tmp26_1 = B_1_2*w36;
2714                                        const double tmp30_1 = B_1_3*w43;
2715                                        const double tmp70_1 = B_0_2*w35;
2716                                        const double tmp34_1 = B_1_3*w39;
2717                                        const double tmp51_1 = B_1_0*w43;
2718                                        const double tmp31_1 = B_0_2*w34;
2719                                        const double tmp45_1 = tmp3_0*w28;
2720                                        const double tmp11_1 = tmp1_0*w39;
2721                                        const double tmp52_1 = B_1_1*w29;
2722                                        const double tmp44_1 = B_1_3*w32;
2723                                        const double tmp25_1 = B_1_1*w39;
2724                                        const double tmp47_1 = tmp2_0*w33;
2725                                        const double tmp72_1 = B_1_3*w29;
2726                                        const double tmp40_1 = tmp2_0*w28;
2727                                        const double tmp46_1 = B_1_2*w40;
2728                                        const double tmp36_1 = B_1_2*w42;
2729                                        const double tmp24_1 = B_0_0*w37;
2730                                        const double tmp77_1 = B_0_3*w35;
2731                                        const double tmp68_1 = B_0_3*w37;
2732                                        const double tmp78_1 = B_0_0*w34;
2733                                        const double tmp12_1 = tmp0_0*w36;
2734                                        const double tmp75_1 = B_1_0*w32;
2735                                        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;
2736                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2737                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2738                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
2739                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2740                                        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;
2741                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
2742                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
2743                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
2744                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2745                                        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;
2746                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2747                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
2748                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
2749                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2750                                        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;
2751                                    }
2752                                }
2753                            } else { // constant data
2754                                for (index_t k=0; k<numEq; k++) {
2755                                    for (index_t m=0; m<numComp; m++) {
2756                                        const double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
2757                                        const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
2758                                        const double tmp6_1 = B_1*w50;
2759                                        const double tmp1_1 = B_1*w45;
2760                                        const double tmp5_1 = B_1*w49;
2761                                        const double tmp4_1 = B_1*w48;
2762                                        const double tmp0_1 = B_0*w44;
2763                                        const double tmp2_1 = B_0*w46;
2764                                        const double tmp7_1 = B_0*w51;
2765                                        const double tmp3_1 = B_0*w47;
2766                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;
2767                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
2768                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;
2769                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
2770                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
2771                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;
2772                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;
2773                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;
2774                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;
2775                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
2776                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;
2777                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
2778                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;
2779                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;
2780                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;
2781                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;
2782                                    }
2783                                }
2784                            }
2785                        }
2786                        ///////////////
2787                        // process C //
2788                        ///////////////
2789                        if (!C.isEmpty()) {
2790                            add_EM_S=true;
2791                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2792                            if (C.actsExpanded()) {
2793                                for (index_t k=0; k<numEq; k++) {
2794                                    for (index_t m=0; m<numComp; m++) {
2795                                        const double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
2796                                        const double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
2797                                        const double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
2798                                        const double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
2799                                        const double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];
2800                                        const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
2801                                        const double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];
2802                                        const double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];
2803                                        const double tmp2_0 = C_0_1 + C_0_3;
2804                                        const double tmp1_0 = C_1_2 + C_1_3;
2805                                        const double tmp3_0 = C_0_0 + C_0_2;
2806                                        const double tmp0_0 = C_1_0 + C_1_1;
2807                                        const double tmp64_1 = C_0_2*w30;
2808                                        const double tmp14_1 = C_0_2*w28;
2809                                        const double tmp19_1 = C_0_3*w31;
2810                                        const double tmp22_1 = C_1_0*w40;
2811                                        const double tmp37_1 = tmp3_0*w35;
2812                                        const double tmp29_1 = C_0_1*w35;
2813                                        const double tmp73_1 = C_0_2*w37;
2814                                        const double tmp74_1 = C_1_2*w41;
2815                                        const double tmp52_1 = C_1_3*w39;
2816                                        const double tmp25_1 = C_1_1*w39;
2817                                        const double tmp62_1 = C_0_1*w31;
2818                                        const double tmp79_1 = C_1_1*w40;
2819                                        const double tmp43_1 = C_1_1*w36;
2820                                        const double tmp27_1 = C_0_3*w38;
2821                                        const double tmp28_1 = C_1_0*w42;
2822                                        const double tmp63_1 = C_1_1*w42;
2823                                        const double tmp59_1 = C_0_3*w34;
2824                                        const double tmp72_1 = C_1_3*w29;
2825                                        const double tmp40_1 = tmp2_0*w35;
2826                                        const double tmp13_1 = C_0_3*w30;
2827                                        const double tmp51_1 = C_1_2*w40;
2828                                        const double tmp54_1 = C_1_2*w42;
2829                                        const double tmp12_1 = C_0_0*w31;
2830                                        const double tmp2_1 = tmp1_0*w32;
2831                                        const double tmp68_1 = C_0_2*w31;
2832                                        const double tmp75_1 = C_1_0*w32;
2833                                        const double tmp49_1 = C_1_1*w41;
2834                                        const double tmp4_1 = C_0_2*w35;
2835                                        const double tmp66_1 = C_0_3*w28;
2836                                        const double tmp56_1 = C_0_1*w37;
2837                                        const double tmp5_1 = C_0_1*w34;
2838                                        const double tmp38_1 = tmp2_0*w34;
2839                                        const double tmp76_1 = C_0_1*w38;
2840                                        const double tmp21_1 = C_0_1*w28;
2841                                        const double tmp69_1 = C_0_1*w30;
2842                                        const double tmp53_1 = C_1_0*w36;
2843                                        const double tmp42_1 = C_1_2*w39;
2844                                        const double tmp32_1 = tmp1_0*w29;
2845                                        const double tmp45_1 = C_1_0*w43;
2846                                        const double tmp33_1 = tmp0_0*w32;
2847                                        const double tmp35_1 = C_1_0*w41;
2848                                        const double tmp26_1 = C_1_2*w36;
2849                                        const double tmp67_1 = C_0_0*w33;
2850                                        const double tmp31_1 = C_0_2*w34;
2851                                        const double tmp20_1 = C_0_0*w30;
2852                                        const double tmp70_1 = C_0_0*w28;
2853                                        const double tmp8_1 = tmp0_0*w39;
2854                                        const double tmp30_1 = C_1_3*w43;
2855                                        const double tmp0_1 = tmp0_0*w29;
2856                                        const double tmp17_1 = C_1_3*w41;
2857                                        const double tmp58_1 = C_0_0*w35;
2858                                        const double tmp9_1 = tmp3_0*w33;
2859                                        const double tmp61_1 = C_1_3*w36;
2860                                        const double tmp41_1 = tmp3_0*w34;
2861                                        const double tmp50_1 = C_1_3*w32;
2862                                        const double tmp18_1 = C_1_1*w32;
2863                                        const double tmp6_1 = tmp1_0*w36;
2864                                        const double tmp3_1 = C_0_0*w38;
2865                                        const double tmp34_1 = C_1_1*w29;
2866                                        const double tmp77_1 = C_0_3*w35;
2867                                        const double tmp65_1 = C_1_2*w43;
2868                                        const double tmp71_1 = C_0_3*w33;
2869                                        const double tmp55_1 = C_1_1*w43;
2870                                        const double tmp46_1 = tmp3_0*w28;
2871                                        const double tmp24_1 = C_0_0*w37;
2872                                        const double tmp10_1 = tmp1_0*w39;
2873                                        const double tmp48_1 = C_1_0*w29;
2874                                        const double tmp15_1 = C_0_1*w33;
2875                                        const double tmp36_1 = C_1_2*w32;
2876                                        const double tmp60_1 = C_1_0*w39;
2877                                        const double tmp47_1 = tmp2_0*w33;
2878                                        const double tmp16_1 = C_1_2*w29;
2879                                        const double tmp1_1 = C_0_3*w37;
2880                                        const double tmp7_1 = tmp2_0*w28;
2881                                        const double tmp39_1 = C_1_3*w40;
2882                                        const double tmp44_1 = C_1_3*w42;
2883                                        const double tmp57_1 = C_0_2*w38;
2884                                        const double tmp78_1 = C_0_0*w34;
2885                                        const double tmp11_1 = tmp0_0*w36;
2886                                        const double tmp23_1 = C_0_2*w33;
2887                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2888                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2889                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2890                                        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;
2891                                        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;
2892                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
2893                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
2894                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
2895                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
2896                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2897                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
2898                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2899                                        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;
2900                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2901                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
2902                                        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;
2903                                    }
2904                                }
2905                            } else { // constant data
2906                                for (index_t k=0; k<numEq; k++) {
2907                                    for (index_t m=0; m<numComp; m++) {
2908                                        const double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
2909                                        const double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
2910                                        const double tmp1_1 = C_1*w45;
2911                                        const double tmp3_1 = C_0*w51;
2912                                        const double tmp4_1 = C_0*w44;
2913                                        const double tmp7_1 = C_0*w46;
2914                                        const double tmp5_1 = C_1*w49;
2915                                        const double tmp2_1 = C_1*w48;
2916                                        const double tmp0_1 = C_0*w47;
2917                                        const double tmp6_1 = C_1*w50;
2918                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
2919                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
2920                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
2921                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp4_1 + tmp5_1;
2922                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp6_1;
2923                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
2924                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp7_1;
2925                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp7_1;
2926                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp3_1 + tmp6_