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