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