/[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 3791 by caltinay, Wed Feb 1 05:10:22 2012 UTC trunk/ripley/src/Rectangle.cpp revision 4622 by sshaw, Fri Jan 17 04:55:41 2014 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2012 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>
 extern "C" {  
17  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
18  }  #include <esysUtils/esysFileWriter.h>
19    #include <ripley/DefaultAssembler2D.h>
20    #include <ripley/WaveAssembler2D.h>
21    #include <boost/scoped_array.hpp>
22    #include "esysUtils/EsysRandom.h"
23    
24    #ifdef USE_NETCDF
25    #include <netcdfcpp.h>
26    #endif
27    
28  #if USE_SILO  #if USE_SILO
29  #include <silo.h>  #include <silo.h>
# Line 26  extern "C" { Line 35  extern "C" {
35  #include <iomanip>  #include <iomanip>
36    
37  using namespace std;  using namespace std;
38    using esysUtils::FileWriter;
39    
40  namespace ripley {  namespace ripley {
41    
42  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
43                       double y1, int d0, int d1) :                       double y1, int d0, int d1,
44      RipleyDomain(2),                       const std::vector<double>& points,
45      m_gNE0(n0),                       const std::vector<int>& tags,
46      m_gNE1(n1),                       const std::map<std::string, int>& tagnamestonums) :
47      m_x0(x0),      RipleyDomain(2)
48      m_y0(y0),  {
49      m_l0(x1-x0),      // ignore subdivision parameters for serial run
50      m_l1(y1-y0),      if (m_mpiInfo->size == 1) {
51      m_NX(d0),          d0=1;
52      m_NY(d1)          d1=1;
53  {      }
54    
55        bool warn=false;
56        // if number of subdivisions is non-positive, try to subdivide by the same
57        // ratio as the number of elements
58        if (d0<=0 && d1<=0) {
59            warn=true;
60            d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
61            d1=m_mpiInfo->size/d0;
62            if (d0*d1 != m_mpiInfo->size) {
63                // ratios not the same so subdivide side with more elements only
64                if (n0>n1) {
65                    d0=0;
66                    d1=1;
67                } else {
68                    d0=1;
69                    d1=0;
70                }
71            }
72        }
73        if (d0<=0) {
74            // d1 is preset, determine d0 - throw further down if result is no good
75            d0=m_mpiInfo->size/d1;
76        } else if (d1<=0) {
77            // d0 is preset, determine d1 - throw further down if result is no good
78            d1=m_mpiInfo->size/d0;
79        }
80    
81      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
82      // among number of ranks      // among number of ranks
83      if (m_NX*m_NY != m_mpiInfo->size)      if (d0*d1 != m_mpiInfo->size)
84          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
85    
86      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)      if (warn) {
87          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
88                << d1 << "). This may not be optimal!" << endl;
89        }
90    
91        double l0 = x1-x0;
92        double l1 = y1-y0;
93        m_dx[0] = l0/n0;
94        m_dx[1] = l1/n1;
95    
96        if ((n0+1)%d0 > 0) {
97            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
98            l0=m_dx[0]*n0;
99            cout << "Warning: Adjusted number of elements and length. N0="
100                << n0 << ", l0=" << l0 << endl;
101        }
102        if ((n1+1)%d1 > 0) {
103            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
104            l1=m_dx[1]*n1;
105            cout << "Warning: Adjusted number of elements and length. N1="
106                << n1 << ", l1=" << l1 << endl;
107        }
108    
109      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))      if ((d0 > 1 && (n0+1)/d0<2) || (d1 > 1 && (n1+1)/d1<2))
110          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
111    
112        m_gNE[0] = n0;
113        m_gNE[1] = n1;
114        m_origin[0] = x0;
115        m_origin[1] = y0;
116        m_length[0] = l0;
117        m_length[1] = l1;
118        m_NX[0] = d0;
119        m_NX[1] = d1;
120    
121      // local number of elements (with and without overlap)      // local number of elements (with and without overlap)
122      m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);      m_NE[0] = m_ownNE[0] = (d0>1 ? (n0+1)/d0 : n0);
123      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)      if (m_mpiInfo->rank%d0>0 && m_mpiInfo->rank%d0<d0-1)
124          m_NE0++;          m_NE[0]++;
125      else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)      else if (d0>1 && m_mpiInfo->rank%d0==d0-1)
126          m_ownNE0--;          m_ownNE[0]--;
127    
128      m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);      m_NE[1] = m_ownNE[1] = (d1>1 ? (n1+1)/d1 : n1);
129      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)      if (m_mpiInfo->rank/d0>0 && m_mpiInfo->rank/d0<d1-1)
130          m_NE1++;          m_NE[1]++;
131      else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)      else if (d1>1 && m_mpiInfo->rank/d0==d1-1)
132          m_ownNE1--;          m_ownNE[1]--;
133    
134      // local number of nodes      // local number of nodes
135      m_N0 = m_NE0+1;      m_NN[0] = m_NE[0]+1;
136      m_N1 = m_NE1+1;      m_NN[1] = m_NE[1]+1;
137    
138      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
139      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);      m_offset[0] = (n0+1)/d0*(m_mpiInfo->rank%d0);
140      if (m_offset0 > 0)      if (m_offset[0] > 0)
141          m_offset0--;          m_offset[0]--;
142      m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);      m_offset[1] = (n1+1)/d1*(m_mpiInfo->rank/d0);
143      if (m_offset1 > 0)      if (m_offset[1] > 0)
144          m_offset1--;          m_offset[1]--;
145    
146      populateSampleIds();      populateSampleIds();
147      createPattern();      createPattern();
148        assembler = new DefaultAssembler2D(this, m_dx, m_NX, m_NE, m_NN);
149        for (map<string, int>::const_iterator i = tagnamestonums.begin();
150                i != tagnamestonums.end(); i++) {
151            setTagMap(i->first, i->second);
152        }
153        addPoints(tags.size(), &points[0], &tags[0]);
154  }  }
155    
156  Rectangle::~Rectangle()  Rectangle::~Rectangle()
157  {  {
158      Paso_SystemMatrixPattern_free(m_pattern);      Paso_SystemMatrixPattern_free(m_pattern);
159      Paso_Connector_free(m_connector);      Paso_Connector_free(m_connector);
160        delete assembler;
161  }  }
162    
163  string Rectangle::getDescription() const  string Rectangle::getDescription() const
# Line 97  bool Rectangle::operator==(const Abstrac Line 170  bool Rectangle::operator==(const Abstrac
170      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
171      if (o) {      if (o) {
172          return (RipleyDomain::operator==(other) &&          return (RipleyDomain::operator==(other) &&
173                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1                  m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1]
174                  && m_x0==o->m_x0 && m_y0==o->m_y0                  && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1]
175                  && m_l0==o->m_l0 && m_l1==o->m_l1                  && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1]
176                  && m_NX==o->m_NX && m_NY==o->m_NY);                  && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1]);
177      }      }
178    
179      return false;      return false;
180  }  }
181    
182    void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
183                const ReaderParameters& params) const
184    {
185    #ifdef USE_NETCDF
186        // check destination function space
187        int myN0, myN1;
188        if (out.getFunctionSpace().getTypeCode() == Nodes) {
189            myN0 = m_NN[0];
190            myN1 = m_NN[1];
191        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
192                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
193            myN0 = m_NE[0];
194            myN1 = m_NE[1];
195        } else
196            throw RipleyException("readNcGrid(): invalid function space for output data object");
197    
198        if (params.first.size() != 2)
199            throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
200    
201        if (params.numValues.size() != 2)
202            throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
203    
204        if (params.multiplier.size() != 2)
205            throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");
206        for (size_t i=0; i<params.multiplier.size(); i++)
207            if (params.multiplier[i]<1)
208                throw RipleyException("readNcGrid(): all multipliers must be positive");
209        if (params.reverse.size() != 2)
210            throw RipleyException("readNcGrid(): argument 'reverse' must have 2 entries");
211    
212        // check file existence and size
213        NcFile f(filename.c_str(), NcFile::ReadOnly);
214        if (!f.is_valid())
215            throw RipleyException("readNcGrid(): cannot open file");
216    
217        NcVar* var = f.get_var(varname.c_str());
218        if (!var)
219            throw RipleyException("readNcGrid(): invalid variable");
220    
221        // TODO: rank>0 data support
222        const int numComp = out.getDataPointSize();
223        if (numComp > 1)
224            throw RipleyException("readNcGrid(): only scalar data supported");
225    
226        const int dims = var->num_dims();
227        boost::scoped_array<long> edges(var->edges());
228    
229        // is this a slice of the data object (dims!=2)?
230        // note the expected ordering of edges (as in numpy: y,x)
231        if ( (dims==2 && (params.numValues[1] > edges[0] || params.numValues[0] > edges[1]))
232                || (dims==1 && params.numValues[1]>1) ) {
233            throw RipleyException("readNcGrid(): not enough data in file");
234        }
235    
236        // check if this rank contributes anything
237        if (params.first[0] >= m_offset[0]+myN0 ||
238                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
239                params.first[1] >= m_offset[1]+myN1 ||
240                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1])
241            return;
242    
243        // now determine how much this rank has to write
244    
245        // first coordinates in data object to write to
246        const int first0 = max(0, params.first[0]-m_offset[0]);
247        const int first1 = max(0, params.first[1]-m_offset[1]);
248        // indices to first value in file (not accounting for reverse yet)
249        int idx0 = max(0, m_offset[0]-params.first[0]);
250        int idx1 = max(0, m_offset[1]-params.first[1]);
251        // number of values to read
252        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
253        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
254    
255        // make sure we read the right block if going backwards through file
256        if (params.reverse[0])
257            idx0 = edges[dims-1]-num0-idx0;
258        if (dims>1 && params.reverse[1])
259            idx1 = edges[dims-2]-num1-idx1;
260    
261        vector<double> values(num0*num1);
262        if (dims==2) {
263            var->set_cur(idx1, idx0);
264            var->get(&values[0], num1, num0);
265        } else {
266            var->set_cur(idx0);
267            var->get(&values[0], num0);
268        }
269    
270        const int dpp = out.getNumDataPointsPerSample();
271        out.requireWrite();
272    
273        // helpers for reversing
274        const int x0 = (params.reverse[0] ? num0-1 : 0);
275        const int x_mult = (params.reverse[0] ? -1 : 1);
276        const int y0 = (params.reverse[1] ? num1-1 : 0);
277        const int y_mult = (params.reverse[1] ? -1 : 1);
278    
279        for (index_t y=0; y<num1; y++) {
280    #pragma omp parallel for
281            for (index_t x=0; x<num0; x++) {
282                const int baseIndex = first0+x*params.multiplier[0]
283                                      +(first1+y*params.multiplier[1])*myN0;
284                const int srcIndex = (y0+y_mult*y)*num0+(x0+x_mult*x);
285                if (!isnan(values[srcIndex])) {
286                    for (index_t m1=0; m1<params.multiplier[1]; m1++) {
287                        for (index_t m0=0; m0<params.multiplier[0]; m0++) {
288                            const int dataIndex = baseIndex+m0+m1*myN0;
289                            double* dest = out.getSampleDataRW(dataIndex);
290                            for (index_t q=0; q<dpp; q++) {
291                                *dest++ = values[srcIndex];
292                            }
293                        }
294                    }
295                }
296            }
297        }
298    #else
299        throw RipleyException("readNcGrid(): not compiled with netCDF support");
300    #endif
301    }
302    
303    void Rectangle::readBinaryGrid(escript::Data& out, string filename,
304                                   const ReaderParameters& params) const
305    {
306        // the mapping is not universally correct but should work on our
307        // supported platforms
308        switch (params.dataType) {
309            case DATATYPE_INT32:
310                readBinaryGridImpl<int>(out, filename, params);
311                break;
312            case DATATYPE_FLOAT32:
313                readBinaryGridImpl<float>(out, filename, params);
314                break;
315            case DATATYPE_FLOAT64:
316                readBinaryGridImpl<double>(out, filename, params);
317                break;
318            default:
319                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
320        }
321    }
322    
323    template<typename ValueType>
324    void Rectangle::readBinaryGridImpl(escript::Data& out, const string& filename,
325                                       const ReaderParameters& params) const
326    {
327        // check destination function space
328        int myN0, myN1;
329        if (out.getFunctionSpace().getTypeCode() == Nodes) {
330            myN0 = m_NN[0];
331            myN1 = m_NN[1];
332        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
333                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
334            myN0 = m_NE[0];
335            myN1 = m_NE[1];
336        } else
337            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
338    
339        // check file existence and size
340        ifstream f(filename.c_str(), ifstream::binary);
341        if (f.fail()) {
342            throw RipleyException("readBinaryGrid(): cannot open file");
343        }
344        f.seekg(0, ios::end);
345        const int numComp = out.getDataPointSize();
346        const int filesize = f.tellg();
347        const int reqsize = params.numValues[0]*params.numValues[1]*numComp*sizeof(ValueType);
348        if (filesize < reqsize) {
349            f.close();
350            throw RipleyException("readBinaryGrid(): not enough data in file");
351        }
352    
353        // check if this rank contributes anything
354        if (params.first[0] >= m_offset[0]+myN0 ||
355                params.first[0]+params.numValues[0] <= m_offset[0] ||
356                params.first[1] >= m_offset[1]+myN1 ||
357                params.first[1]+params.numValues[1] <= m_offset[1]) {
358            f.close();
359            return;
360        }
361    
362        // now determine how much this rank has to write
363    
364        // first coordinates in data object to write to
365        const int first0 = max(0, params.first[0]-m_offset[0]);
366        const int first1 = max(0, params.first[1]-m_offset[1]);
367        // indices to first value in file
368        const int idx0 = max(0, m_offset[0]-params.first[0]);
369        const int idx1 = max(0, m_offset[1]-params.first[1]);
370        // number of values to read
371        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
372        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
373    
374        out.requireWrite();
375        vector<ValueType> values(num0*numComp);
376        const int dpp = out.getNumDataPointsPerSample();
377    
378        for (int y=0; y<num1; y++) {
379            const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]);
380            f.seekg(fileofs*sizeof(ValueType));
381            f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
382            for (int x=0; x<num0; x++) {
383                const int baseIndex = first0+x*params.multiplier[0]
384                                        +(first1+y*params.multiplier[1])*myN0;
385                for (int m1=0; m1<params.multiplier[1]; m1++) {
386                    for (int m0=0; m0<params.multiplier[0]; m0++) {
387                        const int dataIndex = baseIndex+m0+m1*myN0;
388                        double* dest = out.getSampleDataRW(dataIndex);
389                        for (int c=0; c<numComp; c++) {
390                            ValueType val = values[x*numComp+c];
391    
392                            if (params.byteOrder != BYTEORDER_NATIVE) {
393                                char* cval = reinterpret_cast<char*>(&val);
394                                // this will alter val!!
395                                byte_swap32(cval);
396                            }
397                            if (!std::isnan(val)) {
398                                for (int q=0; q<dpp; q++) {
399                                    *dest++ = static_cast<double>(val);
400                                }
401                            }
402                        }
403                    }
404                }
405            }
406        }
407    
408        f.close();
409    }
410    
411    void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,
412                                    int byteOrder, int dataType) const
413    {
414        // the mapping is not universally correct but should work on our
415        // supported platforms
416        switch (dataType) {
417            case DATATYPE_INT32:
418                writeBinaryGridImpl<int>(in, filename, byteOrder);
419                break;
420            case DATATYPE_FLOAT32:
421                writeBinaryGridImpl<float>(in, filename, byteOrder);
422                break;
423            case DATATYPE_FLOAT64:
424                writeBinaryGridImpl<double>(in, filename, byteOrder);
425                break;
426            default:
427                throw RipleyException("writeBinaryGrid(): invalid or unsupported datatype");
428        }
429    }
430    
431    template<typename ValueType>
432    void Rectangle::writeBinaryGridImpl(const escript::Data& in,
433                                        const string& filename, int byteOrder) const
434    {
435        // check function space and determine number of points
436        int myN0, myN1;
437        int totalN0, totalN1;
438        if (in.getFunctionSpace().getTypeCode() == Nodes) {
439            myN0 = m_NN[0];
440            myN1 = m_NN[1];
441            totalN0 = m_gNE[0]+1;
442            totalN1 = m_gNE[1]+1;
443        } else if (in.getFunctionSpace().getTypeCode() == Elements ||
444                    in.getFunctionSpace().getTypeCode() == ReducedElements) {
445            myN0 = m_NE[0];
446            myN1 = m_NE[1];
447            totalN0 = m_gNE[0];
448            totalN1 = m_gNE[1];
449        } else
450            throw RipleyException("writeBinaryGrid(): invalid function space of data object");
451    
452        const int numComp = in.getDataPointSize();
453        const int dpp = in.getNumDataPointsPerSample();
454    
455        if (numComp > 1 || dpp > 1)
456            throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
457    
458        escript::Data* _in = const_cast<escript::Data*>(&in);
459        const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;
460    
461        // from here on we know that each sample consists of one value
462        FileWriter fw;
463        fw.openFile(filename, fileSize);
464        MPIBarrier();
465    
466        for (index_t y=0; y<myN1; y++) {
467            const int fileofs = (m_offset[0]+(m_offset[1]+y)*totalN0)*sizeof(ValueType);
468            ostringstream oss;
469    
470            for (index_t x=0; x<myN0; x++) {
471                const double* sample = _in->getSampleDataRO(y*myN0+x);
472                ValueType fvalue = static_cast<ValueType>(*sample);
473                if (byteOrder == BYTEORDER_NATIVE) {
474                    oss.write((char*)&fvalue, sizeof(fvalue));
475                } else {
476                    char* value = reinterpret_cast<char*>(&fvalue);
477                    oss.write(byte_swap32(value), sizeof(fvalue));
478                }
479            }
480            fw.writeAt(oss, fileofs);
481        }
482        fw.close();
483    }
484    
485  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
486  {  {
487  #if USE_SILO  #if USE_SILO
# Line 114  void Rectangle::dump(const string& fileN Line 490  void Rectangle::dump(const string& fileN
490          fn+=".silo";          fn+=".silo";
491      }      }
492    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
493      int driver=DB_HDF5;          int driver=DB_HDF5;    
494      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
495        const char* blockDirFmt = "/block%04d";
496    
497  #ifdef ESYS_MPI  #ifdef ESYS_MPI
498      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
499        const int NUM_SILO_FILES = 1;
500  #endif  #endif
501    
502      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 168  void Rectangle::dump(const string& fileN Line 544  void Rectangle::dump(const string& fileN
544      }      }
545      */      */
546    
547      boost::scoped_ptr<double> x(new double[m_N0]);      boost::scoped_ptr<double> x(new double[m_NN[0]]);
548      boost::scoped_ptr<double> y(new double[m_N1]);      boost::scoped_ptr<double> y(new double[m_NN[1]]);
549      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);  
550  #pragma omp parallel  #pragma omp parallel
551      {      {
552  #pragma omp for nowait  #pragma omp for nowait
553          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
554              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=getLocalCoordinate(i0, 0);
555          }          }
556  #pragma omp for nowait  #pragma omp for nowait
557          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
558              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=getLocalCoordinate(i1, 1);
559          }          }
560      }      }
561      IndexVector dims = getNumNodesPerDim();      int* dims = const_cast<int*>(getNumNodesPerDim());
562    
563      // write mesh      // write mesh
564      DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,
565              DB_COLLINEAR, NULL);              DB_COLLINEAR, NULL);
566    
567      // write node ids      // write node ids
568      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 2,      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 2,
569              NULL, 0, DB_INT, DB_NODECENT, NULL);              NULL, 0, DB_INT, DB_NODECENT, NULL);
570    
571      // write element ids      // write element ids
572      dims = getNumElementsPerDim();      dims = const_cast<int*>(getNumElementsPerDim());
573      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
574              &dims[0], 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);              dims, 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);
575    
576      // rank 0 writes multimesh and multivar      // rank 0 writes multimesh and multivar
577      if (m_mpiInfo->rank == 0) {      if (m_mpiInfo->rank == 0) {
# Line 289  bool Rectangle::ownSample(int fsType, in Line 663  bool Rectangle::ownSample(int fsType, in
663          case Elements:          case Elements:
664          case ReducedElements:          case ReducedElements:
665              // check ownership of element's bottom left node              // check ownership of element's bottom left node
666              return (m_dofMap[id%m_NE0+m_N0*(id/m_NE0)] < getNumDOF());              return (m_dofMap[id%m_NE[0]+m_NN[0]*(id/m_NE[0])] < getNumDOF());
667          case FaceElements:          case FaceElements:
668          case ReducedFaceElements:          case ReducedFaceElements:
669              {              {
670                  // determine which face the sample belongs to before                  // determine which face the sample belongs to before
671                  // checking ownership of corresponding element's first node                  // checking ownership of corresponding element's first node
                 const IndexVector faces = getNumFacesPerBoundary();  
672                  dim_t n=0;                  dim_t n=0;
673                  for (size_t i=0; i<faces.size(); i++) {                  for (size_t i=0; i<4; i++) {
674                      n+=faces[i];                      n+=m_faceCount[i];
675                      if (id<n) {                      if (id<n) {
676                          index_t k;                          index_t k;
677                          if (i==1)                          if (i==1)
678                              k=m_N0-2;                              k=m_NN[0]-2;
679                          else if (i==3)                          else if (i==3)
680                              k=m_N0*(m_N1-2);                              k=m_NN[0]*(m_NN[1]-2);
681                          else                          else
682                              k=0;                              k=0;
683                          // determine whether to move right or up                          // determine whether to move right or up
684                          const index_t delta=(i/2==0 ? m_N0 : 1);                          const index_t delta=(i/2==0 ? m_NN[0] : 1);
685                          return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());                          return (m_dofMap[k+(id-n+m_faceCount[i])*delta] < getNumDOF());
686                      }                      }
687                  }                  }
688                  return false;                  return false;
# Line 331  void Rectangle::setToNormal(escript::Dat Line 704  void Rectangle::setToNormal(escript::Dat
704          {          {
705              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
706  #pragma omp for nowait  #pragma omp for nowait
707                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
708                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
709                      // set vector at two quadrature points                      // set vector at two quadrature points
710                      *o++ = -1.;                      *o++ = -1.;
# Line 343  void Rectangle::setToNormal(escript::Dat Line 716  void Rectangle::setToNormal(escript::Dat
716    
717              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
718  #pragma omp for nowait  #pragma omp for nowait
719                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
720                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
721                      // set vector at two quadrature points                      // set vector at two quadrature points
722                      *o++ = 1.;                      *o++ = 1.;
# Line 355  void Rectangle::setToNormal(escript::Dat Line 728  void Rectangle::setToNormal(escript::Dat
728    
729              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
730  #pragma omp for nowait  #pragma omp for nowait
731                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
732                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
733                      // set vector at two quadrature points                      // set vector at two quadrature points
734                      *o++ = 0.;                      *o++ = 0.;
# Line 367  void Rectangle::setToNormal(escript::Dat Line 740  void Rectangle::setToNormal(escript::Dat
740    
741              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
742  #pragma omp for nowait  #pragma omp for nowait
743                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
744                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
745                      // set vector at two quadrature points                      // set vector at two quadrature points
746                      *o++ = 0.;                      *o++ = 0.;
# Line 383  void Rectangle::setToNormal(escript::Dat Line 756  void Rectangle::setToNormal(escript::Dat
756          {          {
757              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
758  #pragma omp for nowait  #pragma omp for nowait
759                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
760                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
761                      *o++ = -1.;                      *o++ = -1.;
762                      *o = 0.;                      *o = 0.;
# Line 392  void Rectangle::setToNormal(escript::Dat Line 765  void Rectangle::setToNormal(escript::Dat
765    
766              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
767  #pragma omp for nowait  #pragma omp for nowait
768                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
769                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
770                      *o++ = 1.;                      *o++ = 1.;
771                      *o = 0.;                      *o = 0.;
# Line 401  void Rectangle::setToNormal(escript::Dat Line 774  void Rectangle::setToNormal(escript::Dat
774    
775              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
776  #pragma omp for nowait  #pragma omp for nowait
777                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
778                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
779                      *o++ = 0.;                      *o++ = 0.;
780                      *o = -1.;                      *o = -1.;
# Line 410  void Rectangle::setToNormal(escript::Dat Line 783  void Rectangle::setToNormal(escript::Dat
783    
784              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
785  #pragma omp for nowait  #pragma omp for nowait
786                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
787                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
788                      *o++ = 0.;                      *o++ = 0.;
789                      *o = 1.;                      *o = 1.;
# Line 432  void Rectangle::setToSize(escript::Data& Line 805  void Rectangle::setToSize(escript::Data&
805              || out.getFunctionSpace().getTypeCode() == ReducedElements) {              || out.getFunctionSpace().getTypeCode() == ReducedElements) {
806          out.requireWrite();          out.requireWrite();
807          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
808          const double hSize=getFirstCoordAndSpacing(0).second;          const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]);
         const double vSize=getFirstCoordAndSpacing(1).second;  
         const double size=min(hSize,vSize);  
809  #pragma omp parallel for  #pragma omp parallel for
810          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
811              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 444  void Rectangle::setToSize(escript::Data& Line 815  void Rectangle::setToSize(escript::Data&
815              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
816          out.requireWrite();          out.requireWrite();
817          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
         const double hSize=getFirstCoordAndSpacing(0).second;  
         const double vSize=getFirstCoordAndSpacing(1).second;  
818  #pragma omp parallel  #pragma omp parallel
819          {          {
820              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
821  #pragma omp for nowait  #pragma omp for nowait
822                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
823                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
824                      fill(o, o+numQuad, vSize);                      fill(o, o+numQuad, m_dx[1]);
825                  }                  }
826              }              }
827    
828              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
829  #pragma omp for nowait  #pragma omp for nowait
830                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
831                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
832                      fill(o, o+numQuad, vSize);                      fill(o, o+numQuad, m_dx[1]);
833                  }                  }
834              }              }
835    
836              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
837  #pragma omp for nowait  #pragma omp for nowait
838                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
839                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
840                      fill(o, o+numQuad, hSize);                      fill(o, o+numQuad, m_dx[0]);
841                  }                  }
842              }              }
843    
844              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
845  #pragma omp for nowait  #pragma omp for nowait
846                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
847                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
848                      fill(o, o+numQuad, hSize);                      fill(o, o+numQuad, m_dx[0]);
849                  }                  }
850              }              }
851          } // end of parallel section          } // end of parallel section
# Line 489  void Rectangle::setToSize(escript::Data& Line 858  void Rectangle::setToSize(escript::Data&
858      }      }
859  }  }
860    
 Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  
                                                 bool reducedColOrder) const  
 {  
     /* FIXME: reduced  
     if (reducedRowOrder || reducedColOrder)  
         throw RipleyException("getPattern() not implemented for reduced order");  
     */  
     return m_pattern;  
 }  
   
861  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
862  {  {
863      RipleyDomain::Print_Mesh_Info(full);      RipleyDomain::Print_Mesh_Info(full);
# Line 506  void Rectangle::Print_Mesh_Info(const bo Line 865  void Rectangle::Print_Mesh_Info(const bo
865          cout << "     Id  Coordinates" << endl;          cout << "     Id  Coordinates" << endl;
866          cout.precision(15);          cout.precision(15);
867          cout.setf(ios::scientific, ios::floatfield);          cout.setf(ios::scientific, ios::floatfield);
         pair<double,double> xdx = getFirstCoordAndSpacing(0);  
         pair<double,double> ydy = getFirstCoordAndSpacing(1);  
868          for (index_t i=0; i < getNumNodes(); i++) {          for (index_t i=0; i < getNumNodes(); i++) {
869              cout << "  " << setw(5) << m_nodeId[i]              cout << "  " << setw(5) << m_nodeId[i]
870                  << "  " << xdx.first+(i%m_N0)*xdx.second                  << "  " << getLocalCoordinate(i%m_NN[0], 0)
871                  << "  " << ydy.first+(i/m_N0)*ydy.second << endl;                  << "  " << getLocalCoordinate(i/m_NN[0], 1) << endl;
872          }          }
873      }      }
874  }  }
875    
 IndexVector Rectangle::getNumNodesPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_N0);  
     ret.push_back(m_N1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumElementsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NE0);  
     ret.push_back(m_NE1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumFacesPerBoundary() const  
 {  
     IndexVector ret(4, 0);  
     //left  
     if (m_offset0==0)  
         ret[0]=m_NE1;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         ret[1]=m_NE1;  
     //bottom  
     if (m_offset1==0)  
         ret[2]=m_NE0;  
     //top  
     if (m_mpiInfo->rank/m_NX==m_NY-1)  
         ret[3]=m_NE0;  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumSubdivisionsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NX);  
     ret.push_back(m_NY);  
     return ret;  
 }  
   
 pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  
 {  
     if (dim==0) {  
         return pair<double,double>(m_x0+(m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);  
     } else if (dim==1) {  
         return pair<double,double>(m_y0+(m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);  
     }  
     throw RipleyException("getFirstCoordAndSpacing: invalid argument");  
 }  
   
 //protected  
 dim_t Rectangle::getNumDOF() const  
 {  
     return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;  
 }  
   
 //protected  
 dim_t Rectangle::getNumFaceElements() const  
 {  
     const IndexVector faces = getNumFacesPerBoundary();  
     dim_t n=0;  
     for (size_t i=0; i<faces.size(); i++)  
         n+=faces[i];  
     return n;  
 }  
876    
877  //protected  //protected
878  void Rectangle::assembleCoordinates(escript::Data& arg) const  void Rectangle::assembleCoordinates(escript::Data& arg) const
# Line 594  void Rectangle::assembleCoordinates(escr Line 884  void Rectangle::assembleCoordinates(escr
884      if (!numSamplesEqual(&x, 1, getNumNodes()))      if (!numSamplesEqual(&x, 1, getNumNodes()))
885          throw RipleyException("setToX: Illegal number of samples in Data object");          throw RipleyException("setToX: Illegal number of samples in Data object");
886    
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
887      arg.requireWrite();      arg.requireWrite();
888  #pragma omp parallel for  #pragma omp parallel for
889      for (dim_t i1 = 0; i1 < m_N1; i1++) {      for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
890          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
891              double* point = arg.getSampleDataRW(i0+m_N0*i1);              double* point = arg.getSampleDataRW(i0+m_NN[0]*i1);
892              point[0] = xdx.first+i0*xdx.second;              point[0] = getLocalCoordinate(i0, 0);
893              point[1] = ydy.first+i1*ydy.second;              point[1] = getLocalCoordinate(i1, 1);
894          }          }
895      }      }
896  }  }
# Line 611  void Rectangle::assembleCoordinates(escr Line 899  void Rectangle::assembleCoordinates(escr
899  void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const  void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
900  {  {
901      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
902      const double h0 = m_l0/m_gNE0;      const double cx0 = .21132486540518711775/m_dx[0];
903      const double h1 = m_l1/m_gNE1;      const double cx1 = .78867513459481288225/m_dx[0];
904      const double cx0 = -1./h0;      const double cx2 = 1./m_dx[0];
905      const double cx1 = -.78867513459481288225/h0;      const double cy0 = .21132486540518711775/m_dx[1];
906      const double cx2 = -.5/h0;      const double cy1 = .78867513459481288225/m_dx[1];
907      const double cx3 = -.21132486540518711775/h0;      const double cy2 = 1./m_dx[1];
     const double cx4 = .21132486540518711775/h0;  
     const double cx5 = .5/h0;  
     const double cx6 = .78867513459481288225/h0;  
     const double cx7 = 1./h0;  
     const double cy0 = -1./h1;  
     const double cy1 = -.78867513459481288225/h1;  
     const double cy2 = -.5/h1;  
     const double cy3 = -.21132486540518711775/h1;  
     const double cy4 = .21132486540518711775/h1;  
     const double cy5 = .5/h1;  
     const double cy6 = .78867513459481288225/h1;  
     const double cy7 = 1./h1;  
908    
909      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
910          out.requireWrite();          out.requireWrite();
911  #pragma omp parallel for  #pragma omp parallel
912          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
913              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
914                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
915                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
916                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
917                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
918                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
919                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
920                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
921                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
922                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
923                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
924                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
925                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      for (index_t i=0; i < numComp; ++i) {
926                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
927                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
928                  } /* end of component loop i */                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
929              } /* end of k0 loop */                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
930          } /* end of k1 loop */                          o[INDEX3(i,0,2,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
931                            o[INDEX3(i,1,2,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
932                            o[INDEX3(i,0,3,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
933                            o[INDEX3(i,1,3,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
934                        } // end of component loop i
935                    } // end of k0 loop
936                } // end of k1 loop
937            } // end of parallel section
938      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
939          out.requireWrite();          out.requireWrite();
940  #pragma omp parallel for  #pragma omp parallel
941          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
942              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
943                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
944                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
945                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
946                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
947                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
948                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
949                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
950                      o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
951                  } /* end of component loop i */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
952              } /* end of k0 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
953          } /* end of k1 loop */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
954                        for (index_t i=0; i < numComp; ++i) {
955                            o[INDEX3(i,0,0,numComp,2)] = (f_10[i] + f_11[i] - f_00[i] - f_01[i])*cx2/2;
956                            o[INDEX3(i,1,0,numComp,2)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
957                        } // end of component loop i
958                    } // end of k0 loop
959                } // end of k1 loop
960            } // end of parallel section
961      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
962          out.requireWrite();          out.requireWrite();
963  #pragma omp parallel  #pragma omp parallel
964          {          {
965                vector<double> f_00(numComp);
966                vector<double> f_01(numComp);
967                vector<double> f_10(numComp);
968                vector<double> f_11(numComp);
969              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
970  #pragma omp for nowait  #pragma omp for nowait
971                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
972                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
973                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
974                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
975                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
976                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
977                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
978                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
979                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy2;
980                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
981                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy2;
982                      } /* end of component loop i */                      } // end of component loop i
983                  } /* end of k1 loop */                  } // end of k1 loop
984              } /* end of face 0 */              } // end of face 0
985              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
986  #pragma omp for nowait  #pragma omp for nowait
987                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
988                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
989                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
990                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
991                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
992                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
993                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
994                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
995                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_11[i]-f_10[i])*cy2;
996                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
997                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = (f_11[i]-f_10[i])*cy2;
998                      } /* end of component loop i */                      } // end of component loop i
999                  } /* end of k1 loop */                  } // end of k1 loop
1000              } /* end of face 1 */              } // end of face 1
1001              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1002  #pragma omp for nowait  #pragma omp for nowait
1003                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1004                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1005                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
1006                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1007                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
1008                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1009                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1010                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx2;
1011                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
1012                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx2;
1013                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
1014                      } /* end of component loop i */                      } // end of component loop i
1015                  } /* end of k0 loop */                  } // end of k0 loop
1016              } /* end of face 2 */              } // end of face 2
1017              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1018  #pragma omp for nowait  #pragma omp for nowait
1019                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1020                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1021                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1022                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1023                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1024                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1025                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1026                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_11[i]-f_01[i])*cx2;
1027                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
1028                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = (f_11[i]-f_01[i])*cx2;
1029                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
1030                      } /* end of component loop i */                      } // end of component loop i
1031                  } /* end of k0 loop */                  } // end of k0 loop
1032              } /* end of face 3 */              } // end of face 3
1033          } // end of parallel section          } // end of parallel section
1034    
1035      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1036          out.requireWrite();          out.requireWrite();
1037  #pragma omp parallel  #pragma omp parallel
1038          {          {
1039                vector<double> f_00(numComp);
1040                vector<double> f_01(numComp);
1041                vector<double> f_10(numComp);
1042                vector<double> f_11(numComp);
1043              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1044  #pragma omp for nowait  #pragma omp for nowait
1045                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1046                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1047                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1048                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
1049                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
1050                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1051                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1052                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i] + f_11[i] - f_00[i] - f_01[i])*cx2/2;
1053                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy2;
1054                      } /* end of component loop i */                      } // end of component loop i
1055                  } /* end of k1 loop */                  } // end of k1 loop
1056              } /* end of face 0 */              } // end of face 0
1057              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1058  #pragma omp for nowait  #pragma omp for nowait
1059                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1060                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
1061                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
1062                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1063                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1064                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1065                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1066                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i] + f_11[i] - f_00[i] - f_01[i])*cx2/2;
1067                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_11[i]-f_10[i])*cy2;
1068                      } /* end of component loop i */                      } // end of component loop i
1069                  } /* end of k1 loop */                  } // end of k1 loop
1070              } /* end of face 1 */              } // end of face 1
1071              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1072  #pragma omp for nowait  #pragma omp for nowait
1073                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1074                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1075                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
1076                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1077                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
1078                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1079                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1080                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx2;
1081                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
1082                      } /* end of component loop i */                      } // end of component loop i
1083                  } /* end of k0 loop */                  } // end of k0 loop
1084              } /* end of face 2 */              } // end of face 2
1085              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1086  #pragma omp for nowait  #pragma omp for nowait
1087                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1088                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1089                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1090                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1091                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1092                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1093                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1094                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_11[i]-f_01[i])*cx2;
1095                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
1096                      } /* end of component loop i */                      } // end of component loop i
1097                  } /* end of k0 loop */                  } // end of k0 loop
1098              } /* end of face 3 */              } // end of face 3
1099          } // end of parallel section          } // end of parallel section
1100      }      }
1101  }  }
# Line 805  void Rectangle::assembleGradient(escript Line 1104  void Rectangle::assembleGradient(escript
1104  void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const  void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1105  {  {
1106      const dim_t numComp = arg.getDataPointSize();      const dim_t numComp = arg.getDataPointSize();
1107      const double h0 = m_l0/m_gNE0;      const index_t left = (m_offset[0]==0 ? 0 : 1);
1108      const double h1 = m_l1/m_gNE1;      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1109      const index_t left = (m_offset0==0 ? 0 : 1);      const int fs=arg.getFunctionSpace().getTypeCode();
1110      const index_t bottom = (m_offset1==0 ? 0 : 1);      if (fs == Elements && arg.actsExpanded()) {
     if (arg.getFunctionSpace().getTypeCode() == Elements) {  
         const double w = h0*h1/4.;  
1111  #pragma omp parallel  #pragma omp parallel
1112          {          {
1113              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1114                const double w = m_dx[0]*m_dx[1]/4.;
1115  #pragma omp for nowait  #pragma omp for nowait
1116              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1117                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1118                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1119                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1120                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1121                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1122                          const double f2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
1123                          const double f3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
1124                          int_local[i]+=(f0+f1+f2+f3)*w;                          int_local[i]+=(f0+f1+f2+f3)*w;
1125                      }  /* end of component loop i */                      }  // end of component loop i
1126                  } /* end of k0 loop */                  } // end of k0 loop
1127              } /* end of k1 loop */              } // end of k1 loop
   
1128  #pragma omp critical  #pragma omp critical
1129              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1130                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1131          } // end of parallel section          } // end of parallel section
1132      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
1133          const double w = h0*h1;      } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1134            const double w = m_dx[0]*m_dx[1];
1135  #pragma omp parallel  #pragma omp parallel
1136          {          {
1137              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1138  #pragma omp for nowait  #pragma omp for nowait
1139              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1140                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1141                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1142                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1143                          int_local[i]+=f[i]*w;                          int_local[i]+=f[i]*w;
1144                      }  /* end of component loop i */                      }
1145                  } /* end of k0 loop */                  }
1146              } /* end of k1 loop */              }
   
1147  #pragma omp critical  #pragma omp critical
1148              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1149                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1150          } // end of parallel section          } // end of parallel section
1151      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
1152          const double w0 = h0/2.;      } else if (fs == FaceElements && arg.actsExpanded()) {
         const double w1 = h1/2.;  
1153  #pragma omp parallel  #pragma omp parallel
1154          {          {
1155              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1156                const double w0 = m_dx[0]/2.;
1157                const double w1 = m_dx[1]/2.;
1158              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1159  #pragma omp for nowait  #pragma omp for nowait
1160                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1161                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1162                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1163                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1164                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1165                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1166                      }  /* end of component loop i */                      }  // end of component loop i
1167                  } /* end of k1 loop */                  } // end of k1 loop
1168              }              }
1169    
1170              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1171  #pragma omp for nowait  #pragma omp for nowait
1172                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1173                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1174                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1175                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1176                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1177                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1178                      }  /* end of component loop i */                      }  // end of component loop i
1179                  } /* end of k1 loop */                  } // end of k1 loop
1180              }              }
1181    
1182              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1183  #pragma omp for nowait  #pragma omp for nowait
1184                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1185                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1186                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1187                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1188                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1189                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1190                      }  /* end of component loop i */                      }  // end of component loop i
1191                  } /* end of k0 loop */                  } // end of k0 loop
1192              }              }
1193    
1194              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1195  #pragma omp for nowait  #pragma omp for nowait
1196                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1197                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1198                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1199                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1200                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1201                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1202                      }  /* end of component loop i */                      }  // end of component loop i
1203                  } /* end of k0 loop */                  } // end of k0 loop
1204              }              }
   
1205  #pragma omp critical  #pragma omp critical
1206              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1207                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1208          } // end of parallel section          } // end of parallel section
1209      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
1210        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1211  #pragma omp parallel  #pragma omp parallel
1212          {          {
1213              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1214              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1215  #pragma omp for nowait  #pragma omp for nowait
1216                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1217                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1218                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1219                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*m_dx[1];
1220                      }  /* end of component loop i */                      }
1221                  } /* end of k1 loop */                  }
1222              }              }
1223    
1224              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1225  #pragma omp for nowait  #pragma omp for nowait
1226                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1227                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1228                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1229                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*m_dx[1];
1230                      }  /* end of component loop i */                      }
1231                  } /* end of k1 loop */                  }
1232              }              }
1233    
1234              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1235  #pragma omp for nowait  #pragma omp for nowait
1236                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1237                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1238                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1239                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*m_dx[0];
1240                      }  /* end of component loop i */                      }
1241                  } /* end of k0 loop */                  }
1242              }              }
1243    
1244              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1245  #pragma omp for nowait  #pragma omp for nowait
1246                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1247                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1248                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1249                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*m_dx[0];
1250                      }  /* end of component loop i */                      }
1251                  } /* end of k0 loop */                  }
1252              }              }
1253    
1254  #pragma omp critical  #pragma omp critical
1255              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1256                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1257          } // end of parallel section          } // end of parallel section
1258      }      } // function space selector
1259  }  }
1260    
1261  //protected  //protected
1262  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
1263  {  {
1264      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1265      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1266      const int x=node%nDOF0;      const int x=node%nDOF0;
1267      const int y=node/nDOF0;      const int y=node/nDOF0;
1268      dim_t num=0;      dim_t num=0;
# Line 994  void Rectangle::nodesToDOF(escript::Data Line 1292  void Rectangle::nodesToDOF(escript::Data
1292      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1293      out.requireWrite();      out.requireWrite();
1294    
1295      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1296      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1297      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1298      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1299  #pragma omp parallel for  #pragma omp parallel for
1300      for (index_t i=0; i<nDOF1; i++) {      for (index_t i=0; i<nDOF1; i++) {
1301          for (index_t j=0; j<nDOF0; j++) {          for (index_t j=0; j<nDOF0; j++) {
1302              const index_t n=j+left+(i+bottom)*m_N0;              const index_t n=j+left+(i+bottom)*m_NN[0];
1303              const double* src=in.getSampleDataRO(n);              const double* src=in.getSampleDataRO(n);
1304              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
1305          }          }
# Line 1027  void Rectangle::dofToNodes(escript::Data Line 1325  void Rectangle::dofToNodes(escript::Data
1325                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1326          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1327      }      }
1328        Paso_Coupler_free(coupler);
1329  }  }
1330    
1331  //private  //private
1332  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
1333  {  {
1334      // identifiers are ordered from left to right, bottom to top globablly.      // degrees of freedom are numbered from left to right, bottom to top in
1335        // each rank, continuing on the next rank (ranks also go left-right,
1336        // bottom-top).
1337        // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
1338        // helps when writing out data rank after rank.
1339    
1340      // build node distribution vector first.      // build node distribution vector first.
1341      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
1342        // constant for all ranks in this implementation
1343      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1344      const dim_t numDOF=getNumDOF();      const dim_t numDOF=getNumDOF();
1345      for (dim_t k=1; k<m_mpiInfo->size; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
# Line 1045  void Rectangle::populateSampleIds() Line 1349  void Rectangle::populateSampleIds()
1349      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
1350      m_dofId.resize(numDOF);      m_dofId.resize(numDOF);
1351      m_elementId.resize(getNumElements());      m_elementId.resize(getNumElements());
1352    
1353        // populate face element counts
1354        //left
1355        if (m_offset[0]==0)
1356            m_faceCount[0]=m_NE[1];
1357        else
1358            m_faceCount[0]=0;
1359        //right
1360        if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
1361            m_faceCount[1]=m_NE[1];
1362        else
1363            m_faceCount[1]=0;
1364        //bottom
1365        if (m_offset[1]==0)
1366            m_faceCount[2]=m_NE[0];
1367        else
1368            m_faceCount[2]=0;
1369        //top
1370        if (m_mpiInfo->rank/m_NX[0]==m_NX[1]-1)
1371            m_faceCount[3]=m_NE[0];
1372        else
1373            m_faceCount[3]=0;
1374    
1375      m_faceId.resize(getNumFaceElements());      m_faceId.resize(getNumFaceElements());
1376    
1377        const index_t left = (m_offset[0]==0 ? 0 : 1);
1378        const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1379        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1380        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1381    
1382    #define globalNodeId(x,y) \
1383        ((m_offset[0]+x)/nDOF0)*nDOF0*nDOF1+(m_offset[0]+x)%nDOF0 \
1384        + ((m_offset[1]+y)/nDOF1)*nDOF0*nDOF1*m_NX[0]+((m_offset[1]+y)%nDOF1)*nDOF0
1385    
1386        // set corner id's outside the parallel region
1387        m_nodeId[0] = globalNodeId(0, 0);
1388        m_nodeId[m_NN[0]-1] = globalNodeId(m_NN[0]-1, 0);
1389        m_nodeId[m_NN[0]*(m_NN[1]-1)] = globalNodeId(0, m_NN[1]-1);
1390        m_nodeId[m_NN[0]*m_NN[1]-1] = globalNodeId(m_NN[0]-1,m_NN[1]-1);
1391    #undef globalNodeId
1392    
1393  #pragma omp parallel  #pragma omp parallel
1394      {      {
1395          // nodes          // populate degrees of freedom and own nodes (identical id)
1396  #pragma omp for nowait  #pragma omp for nowait
1397          for (dim_t i1=0; i1<m_N1; i1++) {          for (dim_t i=0; i<nDOF1; i++) {
1398              for (dim_t i0=0; i0<m_N0; i0++) {              for (dim_t j=0; j<nDOF0; j++) {
1399                  m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;                  const index_t nodeIdx=j+left+(i+bottom)*m_NN[0];
1400                    const index_t dofIdx=j+i*nDOF0;
1401                    m_dofId[dofIdx] = m_nodeId[nodeIdx]
1402                        = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
1403              }              }
1404          }          }
1405    
1406          // degrees of freedom          // populate the rest of the nodes (shared with other ranks)
1407            if (m_faceCount[0]==0) { // left column
1408    #pragma omp for nowait
1409                for (dim_t i=0; i<nDOF1; i++) {
1410                    const index_t nodeIdx=(i+bottom)*m_NN[0];
1411                    const index_t dofId=(i+1)*nDOF0-1;
1412                    m_nodeId[nodeIdx]
1413                        = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
1414                }
1415            }
1416            if (m_faceCount[1]==0) { // right column
1417    #pragma omp for nowait
1418                for (dim_t i=0; i<nDOF1; i++) {
1419                    const index_t nodeIdx=(i+bottom+1)*m_NN[0]-1;
1420                    const index_t dofId=i*nDOF0;
1421                    m_nodeId[nodeIdx]
1422                        = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
1423                }
1424            }
1425            if (m_faceCount[2]==0) { // bottom row
1426    #pragma omp for nowait
1427                for (dim_t i=0; i<nDOF0; i++) {
1428                    const index_t nodeIdx=i+left;
1429                    const index_t dofId=nDOF0*(nDOF1-1)+i;
1430                    m_nodeId[nodeIdx]
1431                        = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
1432                }
1433            }
1434            if (m_faceCount[3]==0) { // top row
1435  #pragma omp for nowait  #pragma omp for nowait
1436          for (dim_t k=0; k<numDOF; k++)              for (dim_t i=0; i<nDOF0; i++) {
1437              m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;                  const index_t nodeIdx=m_NN[0]*(m_NN[1]-1)+i+left;
1438                    const index_t dofId=i;
1439                    m_nodeId[nodeIdx]
1440                        = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
1441                }
1442            }
1443    
1444          // elements          // populate element id's
1445  #pragma omp for nowait  #pragma omp for nowait
1446          for (dim_t i1=0; i1<m_NE1; i1++) {          for (dim_t i1=0; i1<m_NE[1]; i1++) {
1447              for (dim_t i0=0; i0<m_NE0; i0++) {              for (dim_t i0=0; i0<m_NE[0]; i0++) {
1448                  m_elementId[i0+i1*m_NE0]=(m_offset1+i1)*m_gNE0+m_offset0+i0;                  m_elementId[i0+i1*m_NE[0]]=(m_offset[1]+i1)*m_gNE[0]+m_offset[0]+i0;
1449              }              }
1450          }          }
1451    
# Line 1083  void Rectangle::populateSampleIds() Line 1462  void Rectangle::populateSampleIds()
1462      updateTagsInUse(Elements);      updateTagsInUse(Elements);
1463    
1464      // generate face offset vector and set face tags      // generate face offset vector and set face tags
     const IndexVector facesPerEdge = getNumFacesPerBoundary();  
1465      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
1466      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
1467      m_faceOffset.assign(facesPerEdge.size(), -1);      m_faceOffset.assign(4, -1);
1468      m_faceTags.clear();      m_faceTags.clear();
1469      index_t offset=0;      index_t offset=0;
1470      for (size_t i=0; i<facesPerEdge.size(); i++) {      for (size_t i=0; i<4; i++) {
1471          if (facesPerEdge[i]>0) {          if (m_faceCount[i]>0) {
1472              m_faceOffset[i]=offset;              m_faceOffset[i]=offset;
1473              offset+=facesPerEdge[i];              offset+=m_faceCount[i];
1474              m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);              m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
1475          }          }
1476      }      }
1477      setTagMap("left", LEFT);      setTagMap("left", LEFT);
# Line 1106  void Rectangle::populateSampleIds() Line 1484  void Rectangle::populateSampleIds()
1484  //private  //private
1485  void Rectangle::createPattern()  void Rectangle::createPattern()
1486  {  {
1487      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1488      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1489      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1490      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1491    
1492      // populate node->DOF mapping with own degrees of freedom.      // populate node->DOF mapping with own degrees of freedom.
1493      // The rest is assigned in the loop further down      // The rest is assigned in the loop further down
# Line 1117  void Rectangle::createPattern() Line 1495  void Rectangle::createPattern()
1495  #pragma omp parallel for  #pragma omp parallel for
1496      for (index_t i=bottom; i<bottom+nDOF1; i++) {      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1497          for (index_t j=left; j<left+nDOF0; j++) {          for (index_t j=left; j<left+nDOF0; j++) {
1498              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;              m_dofMap[i*m_NN[0]+j]=(i-bottom)*nDOF0+j-left;
1499          }          }
1500      }      }
1501    
# Line 1130  void Rectangle::createPattern() Line 1508  void Rectangle::createPattern()
1508      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
1509      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
1510      int numShared=0;      int numShared=0;
1511      const int x=m_mpiInfo->rank%m_NX;      const int x=m_mpiInfo->rank%m_NX[0];
1512      const int y=m_mpiInfo->rank/m_NX;      const int y=m_mpiInfo->rank/m_NX[0];
1513      for (int i1=-1; i1<2; i1++) {      for (int i1=-1; i1<2; i1++) {
1514          for (int i0=-1; i0<2; i0++) {          for (int i0=-1; i0<2; i0++) {
1515              // skip this rank              // skip this rank
# Line 1140  void Rectangle::createPattern() Line 1518  void Rectangle::createPattern()
1518              // location of neighbour rank              // location of neighbour rank
1519              const int nx=x+i0;              const int nx=x+i0;
1520              const int ny=y+i1;              const int ny=y+i1;
1521              if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {              if (nx>=0 && ny>=0 && nx<m_NX[0] && ny<m_NX[1]) {
1522                  neighbour.push_back(ny*m_NX+nx);                  neighbour.push_back(ny*m_NX[0]+nx);
1523                  if (i0==0) {                  if (i0==0) {
1524                      // sharing top or bottom edge                      // sharing top or bottom edge
1525                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1526                      const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);                      const int firstNode=(i1==-1 ? left : m_NN[0]*(m_NN[1]-1)+left);
1527                      offsetInShared.push_back(offsetInShared.back()+nDOF0);                      offsetInShared.push_back(offsetInShared.back()+nDOF0);
1528                      for (dim_t i=0; i<nDOF0; i++, numShared++) {                      for (dim_t i=0; i<nDOF0; i++, numShared++) {
1529                          sendShared.push_back(firstDOF+i);                          sendShared.push_back(firstDOF+i);
# Line 1160  void Rectangle::createPattern() Line 1538  void Rectangle::createPattern()
1538                  } else if (i1==0) {                  } else if (i1==0) {
1539                      // sharing left or right edge                      // sharing left or right edge
1540                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1541                      const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);                      const int firstNode=(i0==-1 ? bottom*m_NN[0] : (bottom+1)*m_NN[0]-1);
1542                      offsetInShared.push_back(offsetInShared.back()+nDOF1);                      offsetInShared.push_back(offsetInShared.back()+nDOF1);
1543                      for (dim_t i=0; i<nDOF1; i++, numShared++) {                      for (dim_t i=0; i<nDOF1; i++, numShared++) {
1544                          sendShared.push_back(firstDOF+i*nDOF0);                          sendShared.push_back(firstDOF+i*nDOF0);
# Line 1170  void Rectangle::createPattern() Line 1548  void Rectangle::createPattern()
1548                          colIndices[firstDOF+i*nDOF0].push_back(numShared);                          colIndices[firstDOF+i*nDOF0].push_back(numShared);
1549                          if (i<nDOF1-1)                          if (i<nDOF1-1)
1550                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1551                          m_dofMap[firstNode+i*m_N0]=numDOF+numShared;                          m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
1552                      }                      }
1553                  } else {                  } else {
1554                      // sharing a node                      // sharing a node
1555                      const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);                      const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1556                      const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);                      const int node=(i0+1)/2*(m_NN[0]-1)+(i1+1)/2*m_NN[0]*(m_NN[1]-1);
1557                      offsetInShared.push_back(offsetInShared.back()+1);                      offsetInShared.push_back(offsetInShared.back()+1);
1558                      sendShared.push_back(dof);                      sendShared.push_back(dof);
1559                      recvShared.push_back(numDOF+numShared);                      recvShared.push_back(numDOF+numShared);
# Line 1285  void Rectangle::addToMatrixAndRHS(Paso_S Line 1663  void Rectangle::addToMatrixAndRHS(Paso_S
1663      IndexVector rowIndex;      IndexVector rowIndex;
1664      rowIndex.push_back(m_dofMap[firstNode]);      rowIndex.push_back(m_dofMap[firstNode]);
1665      rowIndex.push_back(m_dofMap[firstNode+1]);      rowIndex.push_back(m_dofMap[firstNode+1]);
1666      rowIndex.push_back(m_dofMap[firstNode+m_N0]);      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]]);
1667      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]+1]);
1668      if (addF) {      if (addF) {
1669          double *F_p=F.getSampleDataRW(0);          double *F_p=F.getSampleDataRW(0);
1670          for (index_t i=0; i<rowIndex.size(); i++) {          for (index_t i=0; i<rowIndex.size(); i++) {
# Line 1309  void Rectangle::interpolateNodesOnElemen Line 1687  void Rectangle::interpolateNodesOnElemen
1687      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1688      if (reduced) {      if (reduced) {
1689          out.requireWrite();          out.requireWrite();
1690          const double c0 = .25;          const double c0 = 0.25;
1691  #pragma omp parallel for  #pragma omp parallel
1692          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1693              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1694                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_01(numComp);
1695                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1696                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              vector<double> f_11(numComp);
1697                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1698                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1699                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1700                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
1701                  } /* end of component loop i */                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
1702              } /* end of k0 loop */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
1703          } /* end of k1 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
1704                        double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1705                        for (index_t i=0; i < numComp; ++i) {
1706                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1707                        } /* end of component loop i */
1708                    } /* end of k0 loop */
1709                } /* end of k1 loop */
1710            } /* end of parallel section */
1711      } else {      } else {
1712          out.requireWrite();          out.requireWrite();
1713          const double c0 = .16666666666666666667;          const double c0 = 0.16666666666666666667;
1714          const double c1 = .044658198738520451079;          const double c1 = 0.044658198738520451079;
1715          const double c2 = .62200846792814621559;          const double c2 = 0.62200846792814621559;
1716  #pragma omp parallel for  #pragma omp parallel
1717          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1718              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1719                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
1720                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1721                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
1722                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1723                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1724                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1725                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
1726                      o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
1727                      o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
1728                      o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
1729                  } /* end of component loop i */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1730              } /* end of k0 loop */                      for (index_t i=0; i < numComp; ++i) {
1731          } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1732                            o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1733                            o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1734                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1735                        } /* end of component loop i */
1736                    } /* end of k0 loop */
1737                } /* end of k1 loop */
1738            } /* end of parallel section */
1739      }      }
1740  }  }
1741    
# Line 1354  void Rectangle::interpolateNodesOnFaces( Line 1746  void Rectangle::interpolateNodesOnFaces(
1746      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1747      if (reduced) {      if (reduced) {
1748          out.requireWrite();          out.requireWrite();
         const double c0 = .5;  
1749  #pragma omp parallel  #pragma omp parallel
1750          {          {
1751                vector<double> f_00(numComp);
1752                vector<double> f_01(numComp);
1753                vector<double> f_10(numComp);
1754                vector<double> f_11(numComp);
1755              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1756  #pragma omp for nowait  #pragma omp for nowait
1757                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1758                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1759                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1760                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1761                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1762                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = (f_00[i] + f_01[i])/2;
1763                      } /* end of component loop i */                      } /* end of component loop i */
1764                  } /* end of k1 loop */                  } /* end of k1 loop */
1765              } /* end of face 0 */              } /* end of face 0 */
1766              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1767  #pragma omp for nowait  #pragma omp for nowait
1768                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1769                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1770                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1771                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1772                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1773                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = (f_10[i] + f_11[i])/2;
1774                      } /* end of component loop i */                      } /* end of component loop i */
1775                  } /* end of k1 loop */                  } /* end of k1 loop */
1776              } /* end of face 1 */              } /* end of face 1 */
1777              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1778  #pragma omp for nowait  #pragma omp for nowait
1779                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1780                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1781                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1782                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1783                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1784                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = (f_00[i] + f_10[i])/2;
1785                      } /* end of component loop i */                      } /* end of component loop i */
1786                  } /* end of k0 loop */                  } /* end of k0 loop */
1787              } /* end of face 2 */              } /* end of face 2 */
1788              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1789  #pragma omp for nowait  #pragma omp for nowait
1790                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1791                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1792                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1793                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1794                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1795                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = (f_01[i] + f_11[i])/2;
1796                      } /* end of component loop i */                      } /* end of component loop i */
1797                  } /* end of k0 loop */                  } /* end of k0 loop */
1798              } /* end of face 3 */              } /* end of face 3 */
1799          } // end of parallel section          } /* end of parallel section */
1800      } else {      } else {
1801          out.requireWrite();          out.requireWrite();
1802          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1803          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1804  #pragma omp parallel  #pragma omp parallel
1805          {          {
1806                vector<double> f_00(numComp);
1807                vector<double> f_01(numComp);
1808                vector<double> f_10(numComp);
1809                vector<double> f_11(numComp);
1810              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1811  #pragma omp for nowait  #pragma omp for nowait
1812                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1813                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1814                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1815                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1816                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1817                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_01[i] + c1*f_00[i];
1818                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_01[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_00[i] + c1*f_01[i];
1819                      } /* end of component loop i */                      } /* end of component loop i */
1820                  } /* end of k1 loop */                  } /* end of k1 loop */
1821              } /* end of face 0 */              } /* end of face 0 */
1822              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1823  #pragma omp for nowait  #pragma omp for nowait
1824                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1825                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1826                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1827                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1828                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1829                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = c1*f_10[i] + c0*f_11[i];
1830                          o[INDEX2(i,numComp,1)] = f_10[i]*c0 + f_11[i]*c1;                          o[INDEX2(i,numComp,1)] = c1*f_11[i] + c0*f_10[i];
1831                      } /* end of component loop i */                      } /* end of component loop i */
1832                  } /* end of k1 loop */                  } /* end of k1 loop */
1833              } /* end of face 1 */              } /* end of face 1 */
1834              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1835  #pragma omp for nowait  #pragma omp for nowait
1836                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1837                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1838                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1839                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1840                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1841                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_10[i] + c1*f_00[i];
1842                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_10[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_00[i] + c1*f_10[i];
1843                      } /* end of component loop i */                      } /* end of component loop i */
1844                  } /* end of k0 loop */                  } /* end of k0 loop */
1845              } /* end of face 2 */              } /* end of face 2 */
1846              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1847  #pragma omp for nowait  #pragma omp for nowait
1848                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1849                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1850                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1851                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1852                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1853                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_11[i] + c1*f_01[i];
1854                          o[INDEX2(i,numComp,1)] = f_01[i]*c0 + f_11[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_01[i] + c1*f_11[i];
1855                      } /* end of component loop i */                      } /* end of component loop i */
1856                  } /* end of k0 loop */                  } /* end of k0 loop */
1857              } /* end of face 3 */              } /* end of face 3 */
1858          } // end of parallel section          } /* end of parallel section */
1859      }      }
1860  }  }
1861    
 //protected  
 void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,  
         escript::Data& rhs, const escript::Data& A, const escript::Data& B,  
         const escript::Data& C, const escript::Data& D,  
         const escript::Data& X, const escript::Data& Y) const  
 {  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     const double w0 = -0.1555021169820365539*h1/h0;  
     const double w1 = 0.041666666666666666667;  
     const double w2 = -0.15550211698203655390;  
     const double w3 = 0.041666666666666666667*h0/h1;  
     const double w4 = 0.15550211698203655390;  
     const double w5 = -0.041666666666666666667;  
     const double w6 = -0.01116454968463011277*h1/h0;  
     const double w7 = 0.011164549684630112770;  
     const double w8 = -0.011164549684630112770;  
     const double w9 = -0.041666666666666666667*h1/h0;  
     const double w10 = -0.041666666666666666667*h0/h1;  
     const double w11 = 0.1555021169820365539*h1/h0;  
     const double w12 = 0.1555021169820365539*h0/h1;  
     const double w13 = 0.01116454968463011277*h0/h1;  
     const double w14 = 0.01116454968463011277*h1/h0;  
     const double w15 = 0.041666666666666666667*h1/h0;  
     const double w16 = -0.01116454968463011277*h0/h1;  
     const double w17 = -0.1555021169820365539*h0/h1;  
     const double w18 = -0.33333333333333333333*h1/h0;  
     const double w19 = 0.25;  
     const double w20 = -0.25;  
     const double w21 = 0.16666666666666666667*h0/h1;  
     const double w22 = -0.16666666666666666667*h1/h0;  
     const double w23 = -0.16666666666666666667*h0/h1;  
     const double w24 = 0.33333333333333333333*h1/h0;  
     const double w25 = 0.33333333333333333333*h0/h1;  
     const double w26 = 0.16666666666666666667*h1/h0;  
     const double w27 = -0.33333333333333333333*h0/h1;  
     const double w28 = -0.032861463941450536761*h1;  
     const double w29 = -0.032861463941450536761*h0;  
     const double w30 = -0.12264065304058601714*h1;  
     const double w31 = -0.0023593469594139828636*h1;  
     const double w32 = -0.008805202725216129906*h0;  
     const double w33 = -0.008805202725216129906*h1;  
     const double w34 = 0.032861463941450536761*h1;  
     const double w35 = 0.008805202725216129906*h1;  
     const double w36 = 0.008805202725216129906*h0;  
     const double w37 = 0.0023593469594139828636*h1;  
     const double w38 = 0.12264065304058601714*h1;  
     const double w39 = 0.032861463941450536761*h0;  
     const double w40 = -0.12264065304058601714*h0;  
     const double w41 = -0.0023593469594139828636*h0;  
     const double w42 = 0.0023593469594139828636*h0;  
     const double w43 = 0.12264065304058601714*h0;  
     const double w44 = -0.16666666666666666667*h1;  
     const double w45 = -0.083333333333333333333*h0;  
     const double w46 = 0.083333333333333333333*h1;  
     const double w47 = 0.16666666666666666667*h1;  
     const double w48 = 0.083333333333333333333*h0;  
     const double w49 = -0.16666666666666666667*h0;  
     const double w50 = 0.16666666666666666667*h0;  
     const double w51 = -0.083333333333333333333*h1;  
     const double w52 = 0.025917019497006092316*h0*h1;  
     const double w53 = 0.0018607582807716854616*h0*h1;  
     const double w54 = 0.0069444444444444444444*h0*h1;  
     const double w55 = 0.09672363354357992482*h0*h1;  
     const double w56 = 0.00049858867864229740201*h0*h1;  
     const double w57 = 0.055555555555555555556*h0*h1;  
     const double w58 = 0.027777777777777777778*h0*h1;  
     const double w59 = 0.11111111111111111111*h0*h1;  
     const double w60 = -0.19716878364870322056*h1;  
     const double w61 = -0.19716878364870322056*h0;  
     const double w62 = -0.052831216351296779436*h0;  
     const double w63 = -0.052831216351296779436*h1;  
     const double w64 = 0.19716878364870322056*h1;  
     const double w65 = 0.052831216351296779436*h1;  
     const double w66 = 0.19716878364870322056*h0;  
     const double w67 = 0.052831216351296779436*h0;  
     const double w68 = -0.5*h1;  
     const double w69 = -0.5*h0;  
     const double w70 = 0.5*h1;  
     const double w71 = 0.5*h0;  
     const double w72 = 0.1555021169820365539*h0*h1;  
     const double w73 = 0.041666666666666666667*h0*h1;  
     const double w74 = 0.01116454968463011277*h0*h1;  
     const double w75 = 0.25*h0*h1;  
   
     rhs.requireWrite();  
 #pragma omp parallel  
     {  
         for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring  
 #pragma omp for  
             for (index_t k1=k1_0; k1<m_NE1; k1+=2) {  
                 for (index_t k0=0; k0<m_NE0; ++k0)  {  
                     bool add_EM_S=false;  
                     bool add_EM_F=false;  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     const index_t e = k0 + m_NE0*k1;  
                     ///////////////  
                     // process A //  
                     ///////////////  
                     if (!A.isEmpty()) {  
                         add_EM_S=true;  
                         const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);  
                         if (A.actsExpanded()) {  
                             const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];  
                             const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];  
                             const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];  
                             const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];  
                             const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];  
                             const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];  
                             const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];  
                             const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];  
                             const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];  
                             const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];  
                             const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];  
                             const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];  
                             const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];  
                             const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];  
                             const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];  
                             const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];  
                             const double tmp0_0 = A_01_0 + A_01_3;  
                             const double tmp1_0 = A_00_0 + A_00_1;  
                             const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;  
                             const double tmp3_0 = A_00_2 + A_00_3;  
                             const double tmp4_0 = A_10_1 + A_10_2;  
                             const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;  
                             const double tmp6_0 = A_01_3 + A_10_0;  
                             const double tmp7_0 = A_01_0 + A_10_3;  
                             const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;  
                             const double tmp9_0 = A_01_0 + A_10_0;  
                             const double tmp12_0 = A_11_0 + A_11_2;  
                             const double tmp10_0 = A_01_3 + A_10_3;  
                             const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;  
                             const double tmp13_0 = A_01_2 + A_10_1;  
                             const double tmp11_0 = A_11_1 + A_11_3;  
                             const double tmp18_0 = A_01_1 + A_10_1;  
                             const double tmp15_0 = A_01_1 + A_10_2;  
                             const double tmp16_0 = A_10_0 + A_10_3;  
                             const double tmp17_0 = A_01_1 + A_01_2;  
                             const double tmp19_0 = A_01_2 + A_10_2;  
                             const double tmp0_1 = A_10_3*w8;  
                             const double tmp1_1 = tmp0_0*w1;  
                             const double tmp2_1 = A_01_1*w4;  
                             const double tmp3_1 = tmp1_0*w0;  
                             const double tmp4_1 = A_01_2*w7;  
                             const double tmp5_1 = tmp2_0*w3;  
                             const double tmp6_1 = tmp3_0*w6;  
                             const double tmp7_1 = A_10_0*w2;  
                             const double tmp8_1 = tmp4_0*w5;  
                             const double tmp9_1 = tmp2_0*w10;  
                             const double tmp14_1 = A_10_0*w8;  
                             const double tmp23_1 = tmp3_0*w14;  
                             const double tmp35_1 = A_01_0*w8;  
                             const double tmp54_1 = tmp13_0*w8;  
                             const double tmp20_1 = tmp9_0*w4;  
                             const double tmp25_1 = tmp12_0*w12;  
                             const double tmp44_1 = tmp7_0*w7;  
                             const double tmp26_1 = tmp10_0*w4;  
                             const double tmp52_1 = tmp18_0*w8;  
                             const double tmp48_1 = A_10_1*w7;  
                             const double tmp46_1 = A_01_3*w8;  
                             const double tmp50_1 = A_01_0*w2;  
                             const double tmp56_1 = tmp19_0*w8;  
                             const double tmp19_1 = A_10_3*w2;  
                             const double tmp47_1 = A_10_2*w4;  
                             const double tmp16_1 = tmp3_0*w0;  
                             const double tmp18_1 = tmp1_0*w6;  
                             const double tmp31_1 = tmp11_0*w12;  
                             const double tmp55_1 = tmp15_0*w2;  
                             const double tmp39_1 = A_10_2*w7;  
                             const double tmp11_1 = tmp6_0*w7;  
                             const double tmp40_1 = tmp11_0*w17;  
                             const double tmp34_1 = tmp15_0*w8;  
                             const double tmp33_1 = tmp14_0*w5;  
                             const double tmp24_1 = tmp11_0*w13;  
                             const double tmp43_1 = tmp17_0*w5;  
                             const double tmp15_1 = A_01_2*w4;  
                             const double tmp53_1 = tmp19_0*w2;  
                             const double tmp27_1 = tmp3_0*w11;  
                             const double tmp32_1 = tmp13_0*w2;  
                             const double tmp10_1 = tmp5_0*w9;  
                             const double tmp37_1 = A_10_1*w4;  
                             const double tmp38_1 = tmp5_0*w15;  
                             const double tmp17_1 = A_01_1*w7;  
                             const double tmp12_1 = tmp7_0*w4;  
                             const double tmp22_1 = tmp10_0*w7;  
                             const double tmp57_1 = tmp18_0*w2;  
                             const double tmp28_1 = tmp9_0*w7;  
                             const double tmp29_1 = tmp1_0*w14;  
                             const double tmp51_1 = tmp11_0*w16;  
                             const double tmp42_1 = tmp12_0*w16;  
                             const double tmp49_1 = tmp12_0*w17;  
                             const double tmp21_1 = tmp1_0*w11;  
                             const double tmp45_1 = tmp6_0*w4;  
                             const double tmp13_1 = tmp8_0*w1;  
                             const double tmp36_1 = tmp16_0*w1;  
                             const double tmp41_1 = A_01_3*w2;  
                             const double tmp30_1 = tmp12_0*w13;  
                             EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                         } else { // constant data  
                             const double A_00 = A_p[INDEX2(0,0,2)];  
                             const double A_10 = A_p[INDEX2(1,0,2)];  
                             const double A_01 = A_p[INDEX2(0,1,2)];  
                             const double A_11 = A_p[INDEX2(1,1,2)];  
                             const double tmp0_0 = A_01 + A_10;  
                             const double tmp0_1 = A_00*w18;  
                             const double tmp1_1 = A_01*w19;  
                             const double tmp2_1 = A_10*w20;  
                             const double tmp3_1 = A_11*w21;  
                             const double tmp4_1 = A_00*w22;  
                             const double tmp5_1 = tmp0_0*w19;  
                             const double tmp6_1 = A_11*w23;  
                             const double tmp7_1 = A_11*w25;  
                             const double tmp8_1 = A_00*w24;  
                             const double tmp9_1 = tmp0_0*w20;  
                             const double tmp10_1 = A_01*w20;  
                             const double tmp11_1 = A_11*w27;  
                             const double tmp12_1 = A_00*w26;  
                             const double tmp13_1 = A_10*w19;  
                             EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
                         }  
                     }  
                     ///////////////  
                     // process B //  
                     ///////////////  
                     if (!B.isEmpty()) {  
                         add_EM_S=true;  
                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);  
                         if (B.actsExpanded()) {  
                             const double B_0_0 = B_p[INDEX2(0,0,2)];  
                             const double B_1_0 = B_p[INDEX2(1,0,2)];  
                             const double B_0_1 = B_p[INDEX2(0,1,2)];  
                             const double B_1_1 = B_p[INDEX2(1,1,2)];  
                             const double B_0_2 = B_p[INDEX2(0,2,2)];  
                             const double B_1_2 = B_p[INDEX2(1,2,2)];  
                             const double B_0_3 = B_p[INDEX2(0,3,2)];  
                             const double B_1_3 = B_p[INDEX2(1,3,2)];  
                             const double tmp0_0 = B_1_0 + B_1_1;  
                             const double tmp1_0 = B_1_2 + B_1_3;  
                             const double tmp2_0 = B_0_1 + B_0_3;  
                             const double tmp3_0 = B_0_0 + B_0_2;  
                             const double tmp63_1 = B_1_1*w42;  
                             const double tmp79_1 = B_1_1*w40;  
                             const double tmp37_1 = tmp3_0*w35;  
                             const double tmp8_1 = tmp0_0*w32;  
                             const double tmp71_1 = B_0_1*w34;  
                             const double tmp19_1 = B_0_3*w31;  
                             const double tmp15_1 = B_0_3*w34;  
                             const double tmp9_1 = tmp3_0*w34;  
                             const double tmp35_1 = B_1_0*w36;  
                             const double tmp66_1 = B_0_3*w28;  
                             const double tmp28_1 = B_1_0*w42;  
                             const double tmp22_1 = B_1_0*w40;  
                             const double tmp16_1 = B_1_2*w29;  
                             const double tmp6_1 = tmp2_0*w35;  
                             const double tmp55_1 = B_1_3*w40;  
                             const double tmp50_1 = B_1_3*w42;  
                             const double tmp7_1 = tmp1_0*w29;  
                             const double tmp1_1 = tmp1_0*w32;  
                             const double tmp57_1 = B_0_3*w30;  
                             const double tmp18_1 = B_1_1*w32;  
                             const double tmp53_1 = B_1_0*w41;  
                             const double tmp61_1 = B_1_3*w36;  
                             const double tmp27_1 = B_0_3*w38;  
                             const double tmp64_1 = B_0_2*w30;  
                             const double tmp76_1 = B_0_1*w38;  
                             const double tmp39_1 = tmp2_0*w34;  
                             const double tmp62_1 = B_0_1*w31;  
                             const double tmp56_1 = B_0_0*w31;  
                             const double tmp49_1 = B_1_1*w36;  
                             const double tmp2_1 = B_0_2*w31;  
                             const double tmp23_1 = B_0_2*w33;  
                             const double tmp38_1 = B_1_1*w43;  
                             const double tmp74_1 = B_1_2*w41;  
                             const double tmp43_1 = B_1_1*w41;  
                             const double tmp58_1 = B_0_2*w28;  
                             const double tmp67_1 = B_0_0*w33;  
                             const double tmp33_1 = tmp0_0*w39;  
                             const double tmp4_1 = B_0_0*w28;  
                             const double tmp20_1 = B_0_0*w30;  
                             const double tmp13_1 = B_0_2*w38;  
                             const double tmp65_1 = B_1_2*w43;  
                             const double tmp0_1 = tmp0_0*w29;  
                             const double tmp41_1 = tmp3_0*w33;  
                             const double tmp73_1 = B_0_2*w37;  
                             const double tmp69_1 = B_0_0*w38;  
                             const double tmp48_1 = B_1_2*w39;  
                             const double tmp59_1 = B_0_1*w33;  
                             const double tmp17_1 = B_1_3*w41;  
                             const double tmp5_1 = B_0_3*w33;  
                             const double tmp3_1 = B_0_1*w30;  
                             const double tmp21_1 = B_0_1*w28;  
                             const double tmp42_1 = B_1_0*w29;  
                             const double tmp54_1 = B_1_2*w32;  
                             const double tmp60_1 = B_1_0*w39;  
                             const double tmp32_1 = tmp1_0*w36;  
                             const double tmp10_1 = B_0_1*w37;  
                             const double tmp14_1 = B_0_0*w35;  
                             const double tmp29_1 = B_0_1*w35;  
                             const double tmp26_1 = B_1_2*w36;  
                             const double tmp30_1 = B_1_3*w43;  
                             const double tmp70_1 = B_0_2*w35;  
                             const double tmp34_1 = B_1_3*w39;  
                             const double tmp51_1 = B_1_0*w43;  
                             const double tmp31_1 = B_0_2*w34;  
                             const double tmp45_1 = tmp3_0*w28;  
                             const double tmp11_1 = tmp1_0*w39;  
                             const double tmp52_1 = B_1_1*w29;  
                             const double tmp44_1 = B_1_3*w32;  
                             const double tmp25_1 = B_1_1*w39;  
                             const double tmp47_1 = tmp2_0*w33;  
                             const double tmp72_1 = B_1_3*w29;  
                             const double tmp40_1 = tmp2_0*w28;  
                             const double tmp46_1 = B_1_2*w40;  
                             const double tmp36_1 = B_1_2*w42;  
                             const double tmp24_1 = B_0_0*w37;  
                             const double tmp77_1 = B_0_3*w35;  
                             const double tmp68_1 = B_0_3*w37;  
                             const double tmp78_1 = B_0_0*w34;  
                             const double tmp12_1 = tmp0_0*w36;  
                             const double tmp75_1 = B_1_0*w32;  
                             EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                         } else { // constant data  
                             const double B_0 = B_p[0];  
                             const double B_1 = B_p[1];  
                             const double tmp0_1 = B_0*w44;  
                             const double tmp1_1 = B_1*w45;  
                             const double tmp2_1 = B_0*w46;  
                             const double tmp3_1 = B_0*w47;  
                             const double tmp4_1 = B_1*w48;  
                             const double tmp5_1 = B_1*w49;  
                             const double tmp6_1 = B_1*w50;  
                             const double tmp7_1 = B_0*w51;  
                             EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;  
                         }  
                     }  
                     ///////////////  
                     // process C //  
                     ///////////////  
                     if (!C.isEmpty()) {  
                         add_EM_S=true;  
                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);  
                         if (C.actsExpanded()) {  
                             const double C_0_0 = C_p[INDEX2(0,0,2)];  
                             const double C_1_0 = C_p[INDEX2(1,0,2)];  
                             const double C_0_1 = C_p[INDEX2(0,1,2)];  
                             const double C_1_1 = C_p[INDEX2(1,1,2)];  
                             const double C_0_2 = C_p[INDEX2(0,2,2)];  
                             const double C_1_2 = C_p[INDEX2(1,2,2)];  
                             const double C_0_3 = C_p[INDEX2(0,3,2)];  
                             const double C_1_3 = C_p[INDEX2(1,3,2)];  
                             const double tmp0_0 = C_1_0 + C_1_1;  
                             const double tmp1_0 = C_1_2 + C_1_3;  
                             const double tmp2_0 = C_0_1 + C_0_3;  
                             const double tmp3_0 = C_0_0 + C_0_2;  
                             const double tmp64_1 = C_0_2*w30;  
                             const double tmp14_1 = C_0_2*w28;  
                             const double tmp19_1 = C_0_3*w31;  
                             const double tmp22_1 = C_1_0*w40;  
                             const double tmp37_1 = tmp3_0*w35;  
                             const double tmp29_1 = C_0_1*w35;  
                             const double tmp73_1 = C_0_2*w37;  
                             const double tmp74_1 = C_1_2*w41;  
                             const double tmp52_1 = C_1_3*w39;  
                             const double tmp25_1 = C_1_1*w39;  
                             const double tmp62_1 = C_0_1*w31;  
                             const double tmp79_1 = C_1_1*w40;  
                             const double tmp43_1 = C_1_1*w36;  
                             const double tmp27_1 = C_0_3*w38;  
                             const double tmp28_1 = C_1_0*w42;  
                             const double tmp63_1 = C_1_1*w42;  
                             const double tmp59_1 = C_0_3*w34;  
                             const double tmp72_1 = C_1_3*w29;  
                             const double tmp40_1 = tmp2_0*w35;  
                             const double tmp13_1 = C_0_3*w30;  
                             const double tmp51_1 = C_1_2*w40;  
                             const double tmp54_1 = C_1_2*w42;  
                             const double tmp12_1 = C_0_0*w31;  
                             const double tmp2_1 = tmp1_0*w32;  
                             const double tmp68_1 = C_0_2*w31;  
                             const double tmp75_1 = C_1_0*w32;  
                             const double tmp49_1 = C_1_1*w41;  
                             const double tmp4_1 = C_0_2*w35;  
                             const double tmp66_1 = C_0_3*w28;  
                             const double tmp56_1 = C_0_1*w37;  
                             const double tmp5_1 = C_0_1*w34;  
                             const double tmp38_1 = tmp2_0*w34;  
                             const double tmp76_1 = C_0_1*w38;  
                             const double tmp21_1 = C_0_1*w28;  
                             const double tmp69_1 = C_0_1*w30;  
                             const double tmp53_1 = C_1_0*w36;  
                             const double tmp42_1 = C_1_2*w39;  
                             const double tmp32_1 = tmp1_0*w29;  
                             const double tmp45_1 = C_1_0*w43;  
                             const double tmp33_1 = tmp0_0*w32;  
                             const double tmp35_1 = C_1_0*w41;  
                             const double tmp26_1 = C_1_2*w36;  
                             const double tmp67_1 = C_0_0*w33;  
                             const double tmp31_1 = C_0_2*w34;  
                             const double tmp20_1 = C_0_0*w30;  
                             const double tmp70_1 = C_0_0*w28;  
                             const double tmp8_1 = tmp0_0*w39;  
                             const double tmp30_1 = C_1_3*w43;  
                             const double tmp0_1 = tmp0_0*w29;  
                             const double tmp17_1 = C_1_3*w41;  
                             const double tmp58_1 = C_0_0*w35;  
                             const double tmp9_1 = tmp3_0*w33;  
                             const double tmp61_1 = C_1_3*w36;  
                             const double tmp41_1 = tmp3_0*w34;  
                             const double tmp50_1 = C_1_3*w32;  
                             const double tmp18_1 = C_1_1*w32;  
                             const double tmp6_1 = tmp1_0*w36;  
                             const double tmp3_1 = C_0_0*w38;  
                             const double tmp34_1 = C_1_1*w29;  
                             const double tmp77_1 = C_0_3*w35;  
                             const double tmp65_1 = C_1_2*w43;  
                             const double tmp71_1 = C_0_3*w33;  
                             const double tmp55_1 = C_1_1*w43;  
                             const double tmp46_1 = tmp3_0*w28;  
                             const double tmp24_1 = C_0_0*w37;  
                             const double tmp10_1 = tmp1_0*w39;  
                             const double tmp48_1 = C_1_0*w29;  
                             const double tmp15_1 = C_0_1*w33;  
                             const double tmp36_1 = C_1_2*w32;  
                             const double tmp60_1 = C_1_0*w39;  
                             const double tmp47_1 = tmp2_0*w33;  
                             const double tmp16_1 = C_1_2*w29;  
                             const double tmp1_1 = C_0_3*w37;  
                             const double tmp7_1 = tmp2_0*w28;  
                             const double tmp39_1 = C_1_3*w40;  
                             const double tmp44_1 = C_1_3*w42;  
                             const double tmp57_1 = C_0_2*w38;  
                             const double tmp78_1 = C_0_0*w34;  
                             const double tmp11_1 = tmp0_0*w36;  
                             const double tmp23_1 = C_0_2*w33;  
                             EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                         } else { // constant data  
                             const double C_0 = C_p[0];  
                             const double C_1 = C_p[1];  
                             const double tmp0_1 = C_0*w47;  
                             const double tmp1_1 = C_1*w45;  
                             const double tmp2_1 = C_1*w48;  
                             const double tmp3_1 = C_0*w51;  
                             const double tmp4_1 = C_0*w44;  
                             const double tmp5_1 = C_1*w49;  
                             const double tmp6_1 = C_1*w50;  
                             const double tmp7_1 = C_0*w46;  
                             EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;  
                         }  
                     }  
                     ///////////////  
                     // process D //  
                     ///////////////  
                     if (!D.isEmpty()) {  
                         add_EM_S=true;  
                         const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);  
                         if (D.actsExpanded()) {  
                             const double D_0 = D_p[0];  
                             const double D_1 = D_p[1];  
                             const double D_2 = D_p[2];  
                             const double D_3 = D_p[3];  
                             const double tmp0_0 = D_0 + D_1;  
                             const double tmp1_0 = D_2 + D_3;  
                             const double tmp2_0 = D_0 + D_1 + D_2 + D_3;  
                             const double tmp3_0 = D_1 + D_2;  
                             const double tmp4_0 = D_1 + D_3;  
                             const double tmp5_0 = D_0 + D_2;  
                             const double tmp6_0 = D_0 + D_3;  
                             const double tmp0_1 = tmp0_0*w52;  
                             const double tmp1_1 = tmp1_0*w53;  
                             const double tmp2_1 = tmp2_0*w54;  
                             const double tmp3_1 = tmp1_0*w52;  
                             const double tmp4_1 = tmp0_0*w53;  
                             const double tmp5_1 = tmp3_0*w54;  
                             const double tmp6_1 = D_0*w55;  
                             const double tmp7_1 = D_3*w56;  
                             const double tmp8_1 = D_3*w55;  
                             const double tmp9_1 = D_0*w56;  
                             const double tmp10_1 = tmp4_0*w52;  
                             const double tmp11_1 = tmp5_0*w53;  
                             const double tmp12_1 = tmp5_0*w52;  
                             const double tmp13_1 = tmp4_0*w53;  
                             const double tmp14_1 = tmp6_0*w54;  
                             const double tmp15_1 = D_2*w55;  
                             const double tmp16_1 = D_1*w56;  
                             const double tmp17_1 = D_1*w55;  
                             const double tmp18_1 = D_2*w56;  
                             EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp2_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp2_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp2_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;  
                         } else { // constant data  
                             const double tmp0_1 = D_p[0]*w57;  
                             const double tmp1_1 = D_p[0]*w58;  
                             const double tmp2_1 = D_p[0]*w59;  
                             EM_S[INDEX2(0,0,4)]+=tmp2_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp1_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp2_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp0_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp0_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp1_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp2_1;  
                         }  
                     }  
                     ///////////////  
                     // process X //  
                     ///////////////  
                     if (!X.isEmpty()) {  
                         add_EM_F=true;  
                         const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);  
                         if (X.actsExpanded()) {  
                             const double X_0_0 = X_p[INDEX2(0,0,2)];  
                             const double X_1_0 = X_p[INDEX2(1,0,2)];  
                             const double X_0_1 = X_p[INDEX2(0,1,2)];  
                             const double X_1_1 = X_p[INDEX2(1,1,2)];  
                             const double X_0_2 = X_p[INDEX2(0,2,2)];  
                             const double X_1_2 = X_p[INDEX2(1,2,2)];  
                             const double X_0_3 = X_p[INDEX2(0,3,2)];  
                             const double X_1_3 = X_p[INDEX2(1,3,2)];  
                             const double tmp0_0 = X_0_2 + X_0_3;  
                             const double tmp1_0 = X_1_1 + X_1_3;  
                             const double tmp2_0 = X_1_0 + X_1_2;  
                             const double tmp3_0 = X_0_0 + X_0_1;  
                             const double tmp0_1 = tmp0_0*w63;  
                             const double tmp1_1 = tmp1_0*w62;  
                             const double tmp2_1 = tmp2_0*w61;  
                             const double tmp3_1 = tmp3_0*w60;  
                             const double tmp4_1 = tmp0_0*w65;  
                             const double tmp5_1 = tmp3_0*w64;  
                             const double tmp6_1 = tmp2_0*w62;  
                             const double tmp7_1 = tmp1_0*w61;  
                             const double tmp8_1 = tmp2_0*w66;  
                             const double tmp9_1 = tmp3_0*w63;  
                             const double tmp10_1 = tmp0_0*w60;  
                             const double tmp11_1 = tmp1_0*w67;  
                             const double tmp12_1 = tmp1_0*w66;  
                             const double tmp13_1 = tmp3_0*w65;  
                             const double tmp14_1 = tmp0_0*w64;  
                             const double tmp15_1 = tmp2_0*w67;  
                             EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;  
                             EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;  
                             EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                         } else { // constant data  
                             const double tmp0_1 = X_p[1]*w69;  
                             const double tmp1_1 = X_p[0]*w68;  
                             const double tmp2_1 = X_p[0]*w70;  
                             const double tmp3_1 = X_p[1]*w71;  
                             EM_F[0]+=tmp0_1 + tmp1_1;  
                             EM_F[1]+=tmp0_1 + tmp2_1;  
                             EM_F[2]+=tmp1_1 + tmp3_1;  
                             EM_F[3]+=tmp2_1 + tmp3_1;  
                         }  
                     }  
                     ///////////////  
                     // process Y //  
                     ///////////////  
                     if (!Y.isEmpty()) {  
                         add_EM_F=true;  
                         const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);  
                         if (Y.actsExpanded()) {  
                             const double Y_0 = Y_p[0];  
                             const double Y_1 = Y_p[1];  
                             const double Y_2 = Y_p[2];  
                             const double Y_3 = Y_p[3];  
                             const double tmp0_0 = Y_1 + Y_2;  
                             const double tmp1_0 = Y_0 + Y_3;  
                             const double tmp0_1 = Y_0*w72;  
                             const double tmp1_1 = tmp0_0*w73;  
                             const double tmp2_1 = Y_3*w74;  
                             const double tmp3_1 = Y_1*w72;  
                             const double tmp4_1 = tmp1_0*w73;  
                             const double tmp5_1 = Y_2*w74;  
                             const double tmp6_1 = Y_2*w72;  
                             const double tmp7_1 = Y_1*w74;  
                             const double tmp8_1 = Y_3*w72;  
                             const double tmp9_1 = Y_0*w74;  
                             EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;  
                             EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;  
                             EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;  
                         } else { // constant data  
                             const double tmp0_1 = Y_p[0]*w75;  
                             EM_F[0]+=tmp0_1;  
                             EM_F[1]+=tmp0_1;  
                             EM_F[2]+=tmp0_1;  
                             EM_F[3]+=tmp0_1;  
                         }  
                     }  
1862    
                     // add to matrix (if add_EM_S) and RHS (if add_EM_F)  
                     const index_t firstNode=m_N0*k1+k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 } // end k0 loop  
             } // end k1 loop  
         } // end of colouring  
     } // end of parallel region  
 }  
1863    
1864  //protected  namespace
 void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,  
         escript::Data& rhs, const escript::Data& A, const escript::Data& B,  
         const escript::Data& C, const escript::Data& D,  
         const escript::Data& X, const escript::Data& Y) const  
1865  {  {
1866      const double h0 = m_l0/m_gNE0;      // Calculates a guassian blur colvolution matrix for 2D
1867      const double h1 = m_l1/m_gNE1;      double* get2DGauss(unsigned radius, double sigma)
     const double w0 = -.25*h1/h0;  
     const double w1 = .25;  
     const double w2 = -.25;  
     const double w3 = .25*h0/h1;  
     const double w4 = -.25*h0/h1;  
     const double w5 = .25*h1/h0;  
     const double w6 = -.125*h1;  
     const double w7 = -.125*h0;  
     const double w8 = .125*h1;  
     const double w9 = .125*h0;  
     const double w10 = .0625*h0*h1;  
     const double w11 = -.5*h1;  
     const double w12 = -.5*h0;  
     const double w13 = .5*h1;  
     const double w14 = .5*h0;  
     const double w15 = .25*h0*h1;  
   
     rhs.requireWrite();  
 #pragma omp parallel  
1868      {      {
1869          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          double* arr=new double[(radius*2+1)*(radius*2+1)];
1870  #pragma omp for          double common=M_1_PI*0.5*1/(sigma*sigma);
1871              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {      double total=0;
1872                  for (index_t k0=0; k0<m_NE0; ++k0)  {      int r=static_cast<int>(radius);
1873                      bool add_EM_S=false;      for (int y=-r;y<=r;++y)
1874                      bool add_EM_F=false;      {
1875                      vector<double> EM_S(4*4, 0);          for (int x=-r;x<=r;++x)
1876                      vector<double> EM_F(4, 0);          {        
1877                      const index_t e = k0 + m_NE0*k1;              arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));
1878                      ///////////////  // cout << (x+y*(r*2+1)) << " " << arr[(x+r)+(y+r)*(r*2+1)] << endl;
1879                      // process A //              total+=arr[(x+r)+(y+r)*(r*2+1)];
1880                      ///////////////          }
1881                      if (!A.isEmpty()) {      }
1882                          add_EM_S=true;      double invtotal=1/total;
1883                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);  //cout << "Inv total is "    << invtotal << endl;
1884                          const double A_00 = A_p[INDEX2(0,0,2)];      for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)
1885                          const double A_10 = A_p[INDEX2(1,0,2)];      {
1886                          const double A_01 = A_p[INDEX2(0,1,2)];          arr[p]*=invtotal;
1887                          const double A_11 = A_p[INDEX2(1,1,2)];  //cout << p << "->" << arr[p] << endl;      
1888                          const double tmp0_0 = A_01 + A_10;      }
1889                          const double tmp0_1 = A_11*w3;      return arr;
1890                          const double tmp1_1 = A_00*w0;      }
1891                          const double tmp2_1 = A_01*w1;      
1892                          const double tmp3_1 = A_10*w2;      // applies conv to source to get a point.
1893                          const double tmp4_1 = tmp0_0*w1;      // (xp, yp) are the coords in the source matrix not the destination matrix
1894                          const double tmp5_1 = A_11*w4;      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)
1895                          const double tmp6_1 = A_00*w5;      {
1896                          const double tmp7_1 = tmp0_0*w2;          size_t bx=xp-radius, by=yp-radius;
1897                          const double tmp8_1 = A_10*w1;      size_t sbase=bx+by*width;
1898                          const double tmp9_1 = A_01*w2;      double result=0;
1899                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;      for (int y=0;y<2*radius+1;++y)
1900                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;      {    
1901                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;          for (int x=0;x<2*radius+1;++x)
1902                          EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;          {
1903                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;              result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];
1904                          EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;          }
1905                          EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;      }
1906                          EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;          return result;      
1907                          EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;      }
                         EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;  
                         EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;  
                     }  
                     ///////////////  
                     // process B //  
                     ///////////////  
                     if (!B.isEmpty()) {  
                         add_EM_S=true;  
                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);  
                         const double tmp2_1 = B_p[0]*w8;  
                         const double tmp0_1 = B_p[0]*w6;  
                         const double tmp3_1 = B_p[1]*w9;  
                         const double tmp1_1 = B_p[1]*w7;  
                         EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;  
                         EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;  
                         EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;  
                         EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;  
                     }  
                     ///////////////  
                     // process C //  
                     ///////////////  
                     if (!C.isEmpty()) {  
                         add_EM_S=true;  
                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);  
                         const double tmp1_1 = C_p[1]*w7;  
                         const double tmp0_1 = C_p[0]*w8;  
                         const double tmp3_1 = C_p[0]*w6;  
                         const double tmp2_1 = C_p[1]*w9;  
                         EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;  
                         EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;  
                         EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;  
                         EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;  
                     }  
                     ///////////////  
                     // process D //  
                     ///////////////  
                     if (!D.isEmpty()) {  
                         add_EM_S=true;  
                         const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);  
                         const double tmp0_1 = D_p[0]*w10;  
                         EM_S[INDEX2(0,0,4)]+=tmp0_1;  
                         EM_S[INDEX2(1,0,4)]+=tmp0_1;  
                         EM_S[INDEX2(2,0,4)]+=tmp0_1;  
                         EM_S[INDEX2(3,0,4)]+=tmp0_1;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1;  
                         EM_S[INDEX2(1,1,4)]+=tmp0_1;  
                         EM_S[INDEX2(2,1,4)]+=tmp0_1;  
                         EM_S[INDEX2(3,1,4)]+=tmp0_1;  
                         EM_S[INDEX2(0,2,4)]+=tmp0_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp0_1;  
                         EM_S[INDEX2(2,2,4)]+=tmp0_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp0_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp0_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp0_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1;  
                     }  
                     ///////////////  
                     // process X //  
                     ///////////////  
                     if (!X.isEmpty()) {  
                         add_EM_F=true;  
                         const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);  
                         const double tmp0_1 = X_p[0]*w11;  
                         const double tmp2_1 = X_p[0]*w13;  
                         const double tmp1_1 = X_p[1]*w12;  
                         const double tmp3_1 = X_p[1]*w14;  
                         EM_F[0]+=tmp0_1 + tmp1_1;  
                         EM_F[1]+=tmp1_1 + tmp2_1;  
                         EM_F[2]+=tmp0_1 + tmp3_1;  
                         EM_F[3]+=tmp2_1 + tmp3_1;  
                     }  
                     ///////////////  
                     // process Y //  
                     ///////////////  
                     if (!Y.isEmpty()) {  
                         add_EM_F=true;  
                         const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);  
                         const double tmp0_1 = Y_p[0]*w15;  
                         EM_F[0]+=tmp0_1;  
                         EM_F[1]+=tmp0_1;  
                         EM_F[2]+=tmp0_1;  
                         EM_F[3]+=tmp0_1;  
                     }  
   
                     // add to matrix (if add_EM_S) and RHS (if add_EM_F)  
                     const index_t firstNode=m_N0*k1+k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 } // end k0 loop  
             } // end k1 loop  
         } // end of colouring  
     } // end of parallel region  
1908  }  }
1909    
 //protected  
 void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,  
         escript::Data& rhs, const escript::Data& A, const escript::Data& B,  
         const escript::Data& C, const escript::Data& D,  
         const escript::Data& X, const escript::Data& Y) const  
 {  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     dim_t numEq, numComp;  
     if (!mat)  
         numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());  
     else {  
         numEq=mat->logical_row_block_size;  
         numComp=mat->logical_col_block_size;  
     }  
   
     const double w0 = -0.1555021169820365539*h1/h0;  
     const double w1 = 0.041666666666666666667;  
     const double w2 = -0.15550211698203655390;  
     const double w3 = 0.041666666666666666667*h0/h1;  
     const double w4 = 0.15550211698203655390;  
     const double w5 = -0.041666666666666666667;  
     const double w6 = -0.01116454968463011277*h1/h0;  
     const double w7 = 0.011164549684630112770;  
     const double w8 = -0.011164549684630112770;  
     const double w9 = -0.041666666666666666667*h1/h0;  
     const double w10 = -0.041666666666666666667*h0/h1;  
     const double w11 = 0.1555021169820365539*h1/h0;  
     const double w12 = 0.1555021169820365539*h0/h1;  
     const double w13 = 0.01116454968463011277*h0/h1;  
     const double w14 = 0.01116454968463011277*h1/h0;  
     const double w15 = 0.041666666666666666667*h1/h0;  
     const double w16 = -0.01116454968463011277*h0/h1;  
     const double w17 = -0.1555021169820365539*h0/h1;  
     const double w18 = -0.33333333333333333333*h1/h0;  
     const double w19 = 0.25000000000000000000;  
     const double w20 = -0.25000000000000000000;  
     const double w21 = 0.16666666666666666667*h0/h1;  
     const double w22 = -0.16666666666666666667*h1/h0;  
     const double w23 = -0.16666666666666666667*h0/h1;  
     const double w24 = 0.33333333333333333333*h1/h0;  
     const double w25 = 0.33333333333333333333*h0/h1;  
     const double w26 = 0.16666666666666666667*h1/h0;  
     const double w27 = -0.33333333333333333333*h0/h1;  
     const double w28 = -0.032861463941450536761*h1;  
     const double w29 = -0.032861463941450536761*h0;  
     const double w30 = -0.12264065304058601714*h1;  
     const double w31 = -0.0023593469594139828636*h1;  
     const double w32 = -0.008805202725216129906*h0;  
     const double w33 = -0.008805202725216129906*h1;  
     const double w34 = 0.032861463941450536761*h1;  
     const double w35 = 0.008805202725216129906*h1;  
     const double w36 = 0.008805202725216129906*h0;  
     const double w37 = 0.0023593469594139828636*h1;  
     const double w38 = 0.12264065304058601714*h1;  
     const double w39 = 0.032861463941450536761*h0;  
     const double w40 = -0.12264065304058601714*h0;  
     const double w41 = -0.0023593469594139828636*h0;  
     const double w42 = 0.0023593469594139828636*h0;  
     const double w43 = 0.12264065304058601714*h0;  
     const double w44 = -0.16666666666666666667*h1;  
     const double w45 = -0.083333333333333333333*h0;  
     const double w46 = 0.083333333333333333333*h1;  
     const double w47 = 0.16666666666666666667*h1;  
     const double w48 = 0.083333333333333333333*h0;  
     const double w49 = -0.16666666666666666667*h0;  
     const double w50 = 0.16666666666666666667*h0;  
     const double w51 = -0.083333333333333333333*h1;  
     const double w52 = 0.025917019497006092316*h0*h1;  
     const double w53 = 0.0018607582807716854616*h0*h1;  
     const double w54 = 0.0069444444444444444444*h0*h1;  
     const double w55 = 0.09672363354357992482*h0*h1;  
     const double w56 = 0.00049858867864229740201*h0*h1;  
     const double w57 = 0.055555555555555555556*h0*h1;  
     const double w58 = 0.027777777777777777778*h0*h1;  
     const double w59 = 0.11111111111111111111*h0*h1;  
     const double w60 = -0.19716878364870322056*h1;  
     const double w61 = -0.19716878364870322056*h0;  
     const double w62 = -0.052831216351296779436*h0;  
     const double w63 = -0.052831216351296779436*h1;  
     const double w64 = 0.19716878364870322056*h1;  
     const double w65 = 0.052831216351296779436*h1;  
     const double w66 = 0.19716878364870322056*h0;  
     const double w67 = 0.052831216351296779436*h0;  
     const double w68 = -0.5*h1;  
     const double w69 = -0.5*h0;  
     const double w70 = 0.5*h1;  
     const double w71 = 0.5*h0;  
     const double w72 = 0.1555021169820365539*h0*h1;  
     const double w73 = 0.041666666666666666667*h0*h1;  
     const double w74 = 0.01116454968463011277*h0*h1;  
     const double w75 = 0.25*h0*h1;  
1910    
     rhs.requireWrite();  
 #pragma omp parallel  
     {  
         for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring  
 #pragma omp for  
             for (index_t k1=k1_0; k1<m_NE1; k1+=2) {  
                 for (index_t k0=0; k0<m_NE0; ++k0)  {  
                     bool add_EM_S=false;  
                     bool add_EM_F=false;  
                     vector<double> EM_S(4*4*numEq*numComp, 0);  
                     vector<double> EM_F(4*numEq, 0);  
                     const index_t e = k0 + m_NE0*k1;  
                     ///////////////  
                     // process A //  
                     ///////////////  
                     if (!A.isEmpty()) {  
                         add_EM_S=true;  
                         const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);  
                         if (A.actsExpanded()) {  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];  
                                     const double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];  
                                     const double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];  
                                     const double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];  
                                     const double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];  
                                     const double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];  
                                     const double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];  
                                     const double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];  
                                     const double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];  
                                     const double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];  
                                     const double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];  
                                     const double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];  
                                     const double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];  
                                     const double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];  
                                     const double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];  
                                     const double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];  
                                     const double tmp0_0 = A_01_0 + A_01_3;  
                                     const double tmp1_0 = A_00_0 + A_00_1;  
                                     const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;  
                                     const double tmp3_0 = A_00_2 + A_00_3;  
                                     const double tmp4_0 = A_10_1 + A_10_2;  
                                     const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;  
                                     const double tmp6_0 = A_01_3 + A_10_0;  
                                     const double tmp7_0 = A_01_0 + A_10_3;  
                                     const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;  
                                     const double tmp9_0 = A_01_0 + A_10_0;  
                                     const double tmp10_0 = A_01_3 + A_10_3;  
                                     const double tmp11_0 = A_11_1 + A_11_3;  
                                     const double tmp12_0 = A_11_0 + A_11_2;  
                                     const double tmp13_0 = A_01_2 + A_10_1;  
                                     const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;  
                                     const double tmp15_0 = A_01_1 + A_10_2;  
                                     const double tmp16_0 = A_10_0 + A_10_3;  
                                     const double tmp17_0 = A_01_1 + A_01_2;  
                                     const double tmp18_0 = A_01_1 + A_10_1;  
                                     const double tmp19_0 = A_01_2 + A_10_2;  
                                     const double tmp0_1 = A_10_3*w8;  
                                     const double tmp1_1 = tmp0_0*w1;  
                                     const double tmp2_1 = A_01_1*w4;  
                                     const double tmp3_1 = tmp1_0*w0;  
                                     const double tmp4_1 = A_01_2*w7;  
                                     const double tmp5_1 = tmp2_0*w3;  
                                     const double tmp6_1 = tmp3_0*w6;  
                                     const double tmp7_1 = A_10_0*w2;  
                                     const double tmp8_1 = tmp4_0*w5;  
                                     const double tmp9_1 = tmp2_0*w10;  
                                     const double tmp10_1 = tmp5_0*w9;  
                                     const double tmp11_1 = tmp6_0*w7;  
                                     const double tmp12_1 = tmp7_0*w4;  
                                     const double tmp13_1 = tmp8_0*w1;  
                                     const double tmp14_1 = A_10_0*w8;  
                                     const double tmp15_1 = A_01_2*w4;  
                                     const double tmp16_1 = tmp3_0*w0;  
                                     const double tmp17_1 = A_01_1*w7;  
                                     const double tmp18_1 = tmp1_0*w6;  
                                     const double tmp19_1 = A_10_3*w2;  
                                     const double tmp20_1 = tmp9_0*w4;  
                                     const double tmp21_1 = tmp1_0*w11;  
                                     const double tmp22_1 = tmp10_0*w7;  
                                     const double tmp23_1 = tmp3_0*w14;  
                                     const double tmp24_1 = tmp11_0*w13;  
                                     const double tmp25_1 = tmp12_0*w12;  
                                     const double tmp26_1 = tmp10_0*w4;  
                                     const double tmp27_1 = tmp3_0*w11;  
                                     const double tmp28_1 = tmp9_0*w7;  
                                     const double tmp29_1 = tmp1_0*w14;  
                                     const double tmp30_1 = tmp12_0*w13;  
                                     const double tmp31_1 = tmp11_0*w12;  
                                     const double tmp32_1 = tmp13_0*w2;  
                                     const double tmp33_1 = tmp14_0*w5;  
                                     const double tmp34_1 = tmp15_0*w8;  
                                     const double tmp35_1 = A_01_0*w8;  
                                     const double tmp36_1 = tmp16_0*w1;  
                                     const double tmp37_1 = A_10_1*w4;  
                                     const double tmp38_1 = tmp5_0*w15;  
                                     const double tmp39_1 = A_10_2*w7;  
                                     const double tmp40_1 = tmp11_0*w17;  
                                     const double tmp41_1 = A_01_3*w2;  
                                     const double tmp42_1 = tmp12_0*w16;  
                                     const double tmp43_1 = tmp17_0*w5;  
                                     const double tmp44_1 = tmp7_0*w7;  
                                     const double tmp45_1 = tmp6_0*w4;  
                                     const double tmp46_1 = A_01_3*w8;  
                                     const double tmp47_1 = A_10_2*w4;  
                                     const double tmp48_1 = A_10_1*w7;  
                                     const double tmp49_1 = tmp12_0*w17;  
                                     const double tmp50_1 = A_01_0*w2;  
                                     const double tmp51_1 = tmp11_0*w16;  
                                     const double tmp52_1 = tmp18_0*w8;  
                                     const double tmp53_1 = tmp19_0*w2;  
                                     const double tmp54_1 = tmp13_0*w8;  
                                     const double tmp55_1 = tmp15_0*w2;  
                                     const double tmp56_1 = tmp19_0*w8;  
                                     const double tmp57_1 = tmp18_0*w2;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                                 }  
                             }  
                         } else { // constant data  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];  
                                     const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];  
                                     const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];  
                                     const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];  
                                     const double tmp0_0 = A_01 + A_10;  
                                     const double tmp0_1 = A_00*w18;  
                                     const double tmp1_1 = A_01*w19;  
                                     const double tmp2_1 = A_10*w20;  
                                     const double tmp3_1 = A_11*w21;  
                                     const double tmp4_1 = A_00*w22;  
                                     const double tmp5_1 = tmp0_0*w19;  
                                     const double tmp6_1 = A_11*w23;  
                                     const double tmp7_1 = A_11*w25;  
                                     const double tmp8_1 = A_00*w24;  
                                     const double tmp9_1 = tmp0_0*w20;  
                                     const double tmp10_1 = A_01*w20;  
                                     const double tmp11_1 = A_11*w27;  
                                     const double tmp12_1 = A_00*w26;  
                                     const double tmp13_1 = A_10*w19;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process B //  
                     ///////////////  
                     if (!B.isEmpty()) {  
                         add_EM_S=true;  
                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);  
                         if (B.actsExpanded()) {  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];  
                                     const double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];  
                                     const double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];  
                                     const double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];  
                                     const double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];  
                                     const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];  
                                     const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];  
                                     const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];  
                                     const double tmp0_0 = B_1_0 + B_1_1;  
                                     const double tmp1_0 = B_1_2 + B_1_3;  
                                     const double tmp2_0 = B_0_1 + B_0_3;  
                                     const double tmp3_0 = B_0_0 + B_0_2;  
                                     const double tmp63_1 = B_1_1*w42;  
                                     const double tmp79_1 = B_1_1*w40;  
                                     const double tmp37_1 = tmp3_0*w35;  
                                     const double tmp8_1 = tmp0_0*w32;  
                                     const double tmp71_1 = B_0_1*w34;  
                                     const double tmp19_1 = B_0_3*w31;  
                                     const double tmp15_1 = B_0_3*w34;  
                                     const double tmp9_1 = tmp3_0*w34;  
                                     const double tmp35_1 = B_1_0*w36;  
                                     const double tmp66_1 = B_0_3*w28;  
                                     const double tmp28_1 = B_1_0*w42;  
                                     const double tmp22_1 = B_1_0*w40;  
                                     const double tmp16_1 = B_1_2*w29;  
                                     const double tmp6_1 = tmp2_0*w35;  
                                     const double tmp55_1 = B_1_3*w40;  
                                     const double tmp50_1 = B_1_3*w42;  
                                     const double tmp7_1 = tmp1_0*w29;  
                                     const double tmp1_1 = tmp1_0*w32;  
                                     const double tmp57_1 = B_0_3*w30;  
                                     const double tmp18_1 = B_1_1*w32;  
                                     const double tmp53_1 = B_1_0*w41;  
                                     const double tmp61_1 = B_1_3*w36;  
                                     const double tmp27_1 = B_0_3*w38;  
                                     const double tmp64_1 = B_0_2*w30;  
                                     const double tmp76_1 = B_0_1*w38;  
                                     const double tmp39_1 = tmp2_0*w34;  
                                     const double tmp62_1 = B_0_1*w31;  
                                     const double tmp56_1 = B_0_0*w31;  
                                     const double tmp49_1 = B_1_1*w36;  
                                     const double tmp2_1 = B_0_2*w31;  
                                     const double tmp23_1 = B_0_2*w33;  
                                     const double tmp38_1 = B_1_1*w43;  
                                     const double tmp74_1 = B_1_2*w41;  
                                     const double tmp43_1 = B_1_1*w41;  
                                     const double tmp58_1 = B_0_2*w28;  
                                     const double tmp67_1 = B_0_0*w33;  
                                     const double tmp33_1 = tmp0_0*w39;  
                                     const double tmp4_1 = B_0_0*w28;  
                                     const double tmp20_1 = B_0_0*w30;  
                                     const double tmp13_1 = B_0_2*w38;  
                                     const double tmp65_1 = B_1_2*w43;  
                                     const double tmp0_1 = tmp0_0*w29;  
                                     const double tmp41_1 = tmp3_0*w33;  
                                     const double tmp73_1 = B_0_2*w37;  
                                     const double tmp69_1 = B_0_0*w38;  
                                     const double tmp48_1 = B_1_2*w39;  
                                     const double tmp59_1 = B_0_1*w33;  
                                     const double tmp17_1 = B_1_3*w41;  
                                     const double tmp5_1 = B_0_3*w33;  
                                     const double tmp3_1 = B_0_1*w30;  
                                     const double tmp21_1 = B_0_1*w28;  
                                     const double tmp42_1 = B_1_0*w29;  
                                     const double tmp54_1 = B_1_2*w32;  
                                     const double tmp60_1 = B_1_0*w39;  
                                     const double tmp32_1 = tmp1_0*w36;  
                                     const double tmp10_1 = B_0_1*w37;  
                                     const double tmp14_1 = B_0_0*w35;  
                                     const double tmp29_1 = B_0_1*w35;  
                                     const double tmp26_1 = B_1_2*w36;  
                                     const double tmp30_1 = B_1_3*w43;  
                                     const double tmp70_1 = B_0_2*w35;  
                                     const double tmp34_1 = B_1_3*w39;  
                                     const double tmp51_1 = B_1_0*w43;  
                                     const double tmp31_1 = B_0_2*w34;  
                                     const double tmp45_1 = tmp3_0*w28;  
                                     const double tmp11_1 = tmp1_0*w39;  
                                     const double tmp52_1 = B_1_1*w29;  
                                     const double tmp44_1 = B_1_3*w32;  
                                     const double tmp25_1 = B_1_1*w39;  
                                     const double tmp47_1 = tmp2_0*w33;  
                                     const double tmp72_1 = B_1_3*w29;  
                                     const double tmp40_1 = tmp2_0*w28;  
                                     const double tmp46_1 = B_1_2*w40;  
                                     const double tmp36_1 = B_1_2*w42;  
                                     const double tmp24_1 = B_0_0*w37;  
                                     const double tmp77_1 = B_0_3*w35;  
                                     const double tmp68_1 = B_0_3*w37;  
                                     const double tmp78_1 = B_0_0*w34;  
                                     const double tmp12_1 = tmp0_0*w36;  
                                     const double tmp75_1 = B_1_0*w32;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                                 }  
                             }  
                         } else { // constant data  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];  
                                     const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];  
                                     const double tmp6_1 = B_1*w50;  
                                     const double tmp1_1 = B_1*w45;  
                                     const double tmp5_1 = B_1*w49;  
                                     const double tmp4_1 = B_1*w48;  
                                     const double tmp0_1 = B_0*w44;  
                                     const double tmp2_1 = B_0*w46;  
                                     const double tmp7_1 = B_0*w51;  
                                     const double tmp3_1 = B_0*w47;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process C //  
                     ///////////////  
                     if (!C.isEmpty()) {  
                         add_EM_S=true;  
                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);  
                         if (C.actsExpanded()) {  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];  
                                     const double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];  
                                     const double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];  
                                     const double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];  
                                     const double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];  
                                     const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];  
                                     const double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];  
                                     const double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];  
                                     const double tmp2_0 = C_0_1 + C_0_3;  
                                     const double tmp1_0 = C_1_2 + C_1_3;  
                                     const double tmp3_0 = C_0_0 + C_0_2;  
                                     const double tmp0_0 = C_1_0 + C_1_1;  
                                     const double tmp64_1 = C_0_2*w30;  
                                     const double tmp14_1 = C_0_2*w28;  
                                     const double tmp19_1 = C_0_3*w31;  
                                     const double tmp22_1 = C_1_0*w40;  
                                     const double tmp37_1 = tmp3_0*w35;  
                                     const double tmp29_1 = C_0_1*w35;  
                                     const double tmp73_1 = C_0_2*w37;  
                                     const double tmp74_1 = C_1_2*w41;  
                                     const double tmp52_1 = C_1_3*w39;  
                                     const double tmp25_1 = C_1_1*w39;  
                                     const double tmp62_1 = C_0_1*w31;  
                                     const double tmp79_1 = C_1_1*w40;  
                                     const double tmp43_1 = C_1_1*w36;  
                                     const double tmp27_1 = C_0_3*w38;  
                                     const double tmp28_1 = C_1_0*w42;  
                                     const double tmp63_1 = C_1_1*w42;  
                                     const double tmp59_1 = C_0_3*w34;  
                                     const double tmp72_1 = C_1_3*w29;  
                                     const double tmp40_1 = tmp2_0*w35;  
                                     const double tmp13_1 = C_0_3*w30;  
                                     const double tmp51_1 = C_1_2*w40;  
                                     const double tmp54_1 = C_1_2*w42;  
                                     const double tmp12_1 = C_0_0*w31;  
                                     const double tmp2_1 = tmp1_0*w32;  
                                     const double tmp68_1 = C_0_2*w31;  
                                     const double tmp75_1 = C_1_0*w32;  
                                     const double tmp49_1 = C_1_1*w41;  
                                     const double tmp4_1 = C_0_2*w35;  
                                     const double tmp66_1 = C_0_3*w28;  
                                     const double tmp56_1 = C_0_1*w37;  
                                     const double tmp5_1 = C_0_1*w34;  
                                     const double tmp38_1 = tmp2_0*w34;  
                                     const double tmp76_1 = C_0_1*w38;  
                                     const double tmp21_1 = C_0_1*w28;  
                                     const double tmp69_1 = C_0_1*w30;  
                                     const double tmp53_1 = C_1_0*w36;  
                                     const double tmp42_1 = C_1_2*w39;  
                                     const double tmp32_1 = tmp1_0*w29;  
                                     const double tmp45_1 = C_1_0*w43;  
                                     const double tmp33_1 = tmp0_0*w32;  
                                     const double tmp35_1 = C_1_0*w41;  
                                     const double tmp26_1 = C_1_2*w36;  
                                     const double tmp67_1 = C_0_0*w33;  
                                     const double tmp31_1 = C_0_2*w34;  
                                     const double tmp20_1 = C_0_0*w30;  
                                     const double tmp70_1 = C_0_0*w28;  
                                     const double tmp8_1 = tmp0_0*w39;  
                                     const double tmp30_1 = C_1_3*w43;  
                                     const double tmp0_1 = tmp0_0*w29;  
                                     const double tmp17_1 = C_1_3*w41;  
                                     const double tmp58_1 = C_0_0*w35;  
                                     const double tmp9_1 = tmp3_0*w33;  
                                     const double tmp61_1 = C_1_3*w36;  
                                     const double tmp41_1 = tmp3_0*w34;  
                                     const double tmp50_1 = C_1_3*w32;  
                                     const double tmp18_1 = C_1_1*w32;  
                                     const double tmp6_1 = tmp1_0*w36;  
                                     const double tmp3_1 = C_0_0*w38;  
                                     const double tmp34_1 = C_1_1*w29;  
                                     const double tmp77_1 = C_0_3*w35;  
                                     const double tmp65_1 = C_1_2*w43;  
                                     const double tmp71_1 = C_0_3*w33;  
                                     const double tmp55_1 = C_1_1*w43;  
                                     const double tmp46_1 = tmp3_0*w28;  
                                     const double tmp24_1 = C_0_0*w37;  
                                     const double tmp10_1 = tmp1_0*w39;  
                                     const double tmp48_1 = C_1_0*w29;  
                                     const double tmp15_1 = C_0_1*w33;  
                                     const double tmp36_1 = C_1_2*w32;  
                                     const double tmp60_1 = C_1_0*w39;  
                                     const double tmp47_1 = tmp2_0*w33;  
                                     const double tmp16_1 = C_1_2*w29;  
                                     const double tmp1_1 = C_0_3*w37;  
                                     const double tmp7_1 = tmp2_0*w28;  
                                     const double tmp39_1 = C_1_3*w40;  
                                     const double tmp44_1 = C_1_3*w42;  
                                     const double tmp57_1 = C_0_2*w38;  
                                     const double tmp78_1 = C_0_0*w34;  
                                     const double tmp11_1 = tmp0_0*w36;  
                                     const double tmp23_1 = C_0_2*w33;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                                 }  
                             }  
                         } else { // constant data  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];  
                                     const double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];  
                                     const double tmp1_1 = C_1*w45;  
                                     const double tmp3_1 = C_0*w51;  
                                     const double tmp4_1 = C_0*w44;  
                                     const double tmp7_1 = C_0*w46;  
                                     const double tmp5_1 = C_1*w49;  
                                     const double tmp2_1 = C_1*w48;  
                                     const double tmp0_1 = C_0*w47;  
                                     const double tmp6_1 = C_1*w50;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp3_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp3_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp6_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp4_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp5_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process D //  
                     ///////////////  
                     if (!D.isEmpty()) {  
                         add_EM_S=true;  
                         const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);  
                         if (D.actsExpanded()) {  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const double D_0 = D_p[INDEX3(k, m, 0, numEq, numComp)];  
                                     const double D_1 = D_p[INDEX3(k, m, 1, numEq, numComp)];  
                                     const double D_2 = D_p[INDEX3(k, m, 2, numEq, numComp)];  
                                     const double D_3 = D_p[INDEX3(k, m, 3, numEq, numComp)];  
                                     const double tmp4_0 = D_1 + D_3;  
                                     const double tmp2_0 = D_0 + D_1 + D_2 + D_3;  
                                     const double tmp5_0 = D_0 + D_2;  
                                     const double tmp0_0 = D_0 + D_1;  
                                     const double tmp6_0 = D_0 + D_3;  
                                     const double tmp1_0 = D_2 + D_3;  
                                     const double tmp3_0 = D_1 + D_2;  
                                     const double tmp16_1 = D_1*w56;  
                                     const double tmp14_1 = tmp6_0*w54;  
                                     const double tmp8_1 = D_3*w55;  
                                     const double tmp2_1 = tmp2_0*w54;  
                                     const double tmp12_1 = tmp5_0*w52;  
                                     const double tmp4_1 = tmp0_0*w53;  
                                     const double tmp3_1 = tmp1_0*w52;  
                                     const double tmp13_1 = tmp4_0*w53;  
                                     const double tmp10_1 = tmp4_0*w52;  
                                     const double tmp1_1 = tmp1_0*w53;  
                                     const double tmp7_1 = D_3*w56;  
                                     const double tmp0_1 = tmp0_0*w52;  
                                     const double tmp11_1 = tmp5_0*w53;  
                                     const double tmp9_1 = D_0*w56;  
                                     const double tmp5_1 = tmp3_0*w54;  
                                     const double tmp18_1 = D_2*w56;  
                                     const double tmp17_1 = D_1*w55;  
                                     const double tmp6_1 = D_0*w55;  
                                     const double tmp15_1 = D_2*w55;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp10_1 + tmp11_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp14_1 + tmp17_1 + tmp18_1;  
                                 }  
                              }  
                         } else { // constant data  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const double D_0 = D_p[INDEX2(k, m, numEq)];  
                                     const double tmp2_1 = D_0*w59;  
                                     const double tmp1_1 = D_0*w58;  
                                     const double tmp0_1 = D_0*w57;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp2_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process X //  
                     ///////////////  
                     if (!X.isEmpty()) {  
                         add_EM_F=true;  
                         const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);  
                         if (X.actsExpanded()) {  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double X_0_0 = X_p[INDEX3(k, 0, 0, numEq, 2)];  
                                 const double X_1_0 = X_p[INDEX3(k, 1, 0, numEq, 2)];  
                                 const double X_0_1 = X_p[INDEX3(k, 0, 1, numEq, 2)];  
                                 const double X_1_1 = X_p[INDEX3(k, 1, 1, numEq, 2)];  
                                 const double X_0_2 = X_p[INDEX3(k, 0, 2, numEq, 2)];  
                                 const double X_1_2 = X_p[INDEX3(k, 1, 2, numEq, 2)];  
                                 const double X_0_3 = X_p[INDEX3(k, 0, 3, numEq, 2)];  
                                 const double X_1_3 = X_p[INDEX3(k, 1, 3, numEq, 2)];  
                                 const double tmp1_0 = X_1_1 + X_1_3;  
                                 const double tmp3_0 = X_0_0 + X_0_1;  
                                 const double tmp2_0 = X_1_0 + X_1_2;  
                                 const double tmp0_0 = X_0_2 + X_0_3;  
                                 const double tmp8_1 = tmp2_0*w66;  
                                 const double tmp5_1 = tmp3_0*w64;  
                                 const double tmp14_1 = tmp0_0*w64;  
                                 const double tmp3_1 = tmp3_0*w60;  
                                 const double tmp9_1 = tmp3_0*w63;  
                                 const double tmp13_1 = tmp3_0*w65;  
                                 const double tmp12_1 = tmp1_0*w66;  
                                 const double tmp10_1 = tmp0_0*w60;  
                                 const double tmp2_1 = tmp2_0*w61;  
                                 const double tmp6_1 = tmp2_0*w62;  
                                 const double tmp4_1 = tmp0_0*w65;  
                                 const double tmp11_1 = tmp1_0*w67;  
                                 const double tmp1_1 = tmp1_0*w62;  
                                 const double tmp7_1 = tmp1_0*w61;  
                                 const double tmp0_1 = tmp0_0*w63;  
                                 const double tmp15_1 = tmp2_0*w67;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                             }  
                         } else { // constant data  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double X_0 = X_p[INDEX2(k, 0, numEq)];  
                                 const double X_1 = X_p[INDEX2(k, 1, numEq)];  
                                 const double tmp0_1 = X_1*w69;  
                                 const double tmp1_1 = X_0*w68;  
                                 const double tmp2_1 = X_0*w70;  
                                 const double tmp3_1 = X_1*w71;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process Y //  
                     ///////////////  
                     if (!Y.isEmpty()) {  
                         add_EM_F=true;  
                         const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);  
                         if (Y.actsExpanded()) {  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double Y_0 = Y_p[INDEX2(k, 0, numEq)];  
                                 const double Y_1 = Y_p[INDEX2(k, 1, numEq)];  
                                 const double Y_2 = Y_p[INDEX2(k, 2, numEq)];  
                                 const double Y_3 = Y_p[INDEX2(k, 3, numEq)];  
                                 const double tmp0_0 = Y_1 + Y_2;  
                                 const double tmp1_0 = Y_0 + Y_3;  
                                 const double tmp0_1 = Y_0*w72;  
                                 const double tmp1_1 = tmp0_0*w73;  
                                 const double tmp2_1 = Y_3*w74;  
                                 const double tmp3_1 = Y_1*w72;  
                                 const double tmp4_1 = tmp1_0*w73;  
                                 const double tmp5_1 = Y_2*w74;  
                                 const double tmp6_1 = Y_2*w72;  
                                 const double tmp7_1 = Y_1*w74;  
                                 const double tmp8_1 = Y_3*w72;  
                                 const double tmp9_1 = Y_0*w74;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;  
                             }  
                         } else { // constant data  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double tmp0_1 = Y_p[k]*w75;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
                             }  
                         }  
                     }  
1911    
1912                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)  /* This routine produces a Data object filled with smoothed random data.
1913                      const index_t firstNode=m_N0*k1+k0;  The dimensions of the rectangle being filled are internal[0] x internal[1] points.
1914                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  A parameter radius  dives the size of the stencil used for the smoothing.
1915                              firstNode, numEq, numComp);  A point on the left hand edge for example, will still require `radius` extra points to the left
1916                  } // end k0 loop  in order to complete the stencil.
             } // end k1 loop  
         } // end of colouring  
     } // end of parallel region  
 }  
1917    
1918  //protected  All local calculation is done on an array called `src`, with
1919  void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,  dimensions = ext[0] * ext[1].
1920          escript::Data& rhs, const escript::Data& A, const escript::Data& B,  Where ext[i]= internal[i]+2*radius.
         const escript::Data& C, const escript::Data& D,  
         const escript::Data& X, const escript::Data& Y) const  
 {  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     dim_t numEq, numComp;  
     if (!mat)  
         numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());  
     else {  
         numEq=mat->logical_row_block_size;  
         numComp=mat->logical_col_block_size;  
     }  
   
     const double w0 = -.25*h1/h0;  
     const double w1 = .25;  
     const double w2 = -.25;  
     const double w3 = .25*h0/h1;  
     const double w4 = -.25*h0/h1;  
     const double w5 = .25*h1/h0;  
     const double w6 = -.125*h1;  
     const double w7 = -.125*h0;  
     const double w8 = .125*h1;  
     const double w9 = .125*h0;  
     const double w10 = .0625*h0*h1;  
     const double w11 = -.5*h1;  
     const double w12 = -.5*h0;  
     const double w13 = .5*h1;  
     const double w14 = .5*h0;  
     const double w15 = .25*h0*h1;  
1921    
1922      rhs.requireWrite();  Now for MPI there is overlap to deal with. We need to share both the overlapping
1923  #pragma omp parallel  values themselves but also the external region.
     {  
         for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring  
 #pragma omp for  
             for (index_t k1=k1_0; k1<m_NE1; k1+=2) {  
                 for (index_t k0=0; k0<m_NE0; ++k0)  {  
                     bool add_EM_S=false;  
                     bool add_EM_F=false;  
                     vector<double> EM_S(4*4*numEq*numComp, 0);  
                     vector<double> EM_F(4*numEq, 0);  
                     const index_t e = k0 + m_NE0*k1;  
                     ///////////////  
                     // process A //  
                     ///////////////  
                     if (!A.isEmpty()) {  
                         add_EM_S=true;  
                         const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             for (index_t m=0; m<numComp; m++) {  
                                 const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];  
                                 const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];  
                                 const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];  
                                 const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];  
                                 const double tmp0_0 = A_01 + A_10;  
                                 const double tmp0_1 = A_11*w3;  
                                 const double tmp1_1 = A_00*w0;  
                                 const double tmp2_1 = A_01*w1;  
                                 const double tmp3_1 = A_10*w2;  
                                 const double tmp4_1 = tmp0_0*w1;  
                                 const double tmp5_1 = A_11*w4;  
                                 const double tmp6_1 = A_00*w5;  
                                 const double tmp7_1 = tmp0_0*w2;  
                                 const double tmp8_1 = A_10*w1;  
                                 const double tmp9_1 = A_01*w2;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;  
                                 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
                                 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                                 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;  
                                 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;  
                                 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;  
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;  
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process B //  
                     ///////////////  
                     if (!B.isEmpty()) {  
                         add_EM_S=true;  
                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             for (index_t m=0; m<numComp; m++) {  
                                 const double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];  
                                 const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];  
                                 const double tmp0_1 = B_0*w6;  
                                 const double tmp1_1 = B_1*w7;  
                                 const double tmp2_1 = B_0*w8;  
                                 const double tmp3_1 = B_1*w9;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process C //  
                     ///////////////  
                     if (!C.isEmpty()) {  
                         add_EM_S=true;  
                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             for (index_t m=0; m<numComp; m++) {  
                                 const double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];  
                                 const double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];  
                                 const double tmp0_1 = C_0*w8;  
                                 const double tmp1_1 = C_1*w7;  
                                 const double tmp2_1 = C_1*w9;  
                                 const double tmp3_1 = C_0*w6;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process D //  
                     ///////////////  
                     if (!D.isEmpty()) {  
                         add_EM_S=true;  
                         const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             for (index_t m=0; m<numComp; m++) {  
                                 const double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process X //  
                     ///////////////  
                     if (!X.isEmpty()) {  
                         add_EM_F=true;  
                         const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             const double X_0 = X_p[INDEX2(k, 0, numEq)];  
                             const double X_1 = X_p[INDEX2(k, 1, numEq)];  
                             const double tmp0_1 = X_0*w11;  
                             const double tmp1_1 = X_1*w12;  
                             const double tmp2_1 = X_0*w13;  
                             const double tmp3_1 = X_1*w14;  
                             EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;  
                             EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;  
                             EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp3_1;  
                             EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;  
                         }  
                     }  
                     ///////////////  
                     // process Y //  
                     ///////////////  
                     if (!Y.isEmpty()) {  
                         add_EM_F=true;  
                         const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             const double tmp0_1 = Y_p[k]*w15;  
                             EM_F[INDEX2(k,0,numEq)]+=tmp0_1;  
                             EM_F[INDEX2(k,1,numEq)]+=tmp0_1;  
                             EM_F[INDEX2(k,2,numEq)]+=tmp0_1;  
                             EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
                         }  
                     }  
1924    
1925                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)  In a hypothetical 1-D case:
                     const index_t firstNode=m_N0*k1+k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
                             firstNode, numEq, numComp);  
                 } // end k0 loop  
             } // end k1 loop  
         } // end of colouring  
     } // end of parallel region  
 }  
1926    
 //protected  
 void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,  
       escript::Data& rhs, const escript::Data& d, const escript::Data& y) const  
 {  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     const double w0 = 0.31100423396407310779*h1;  
     const double w1 = 0.022329099369260225539*h1;  
     const double w10 = 0.022329099369260225539*h0;  
     const double w11 = 0.16666666666666666667*h0;  
     const double w12 = 0.33333333333333333333*h0;  
     const double w13 = 0.39433756729740644113*h0;  
     const double w14 = 0.10566243270259355887*h0;  
     const double w15 = 0.5*h0;  
     const double w2 = 0.083333333333333333333*h1;  
     const double w3 = 0.33333333333333333333*h1;  
     const double w4 = 0.16666666666666666667*h1;  
     const double w5 = 0.39433756729740644113*h1;  
     const double w6 = 0.10566243270259355887*h1;  
     const double w7 = 0.5*h1;  
     const double w8 = 0.083333333333333333333*h0;  
     const double w9 = 0.31100423396407310779*h0;  
     const bool add_EM_S=!d.isEmpty();  
     const bool add_EM_F=!y.isEmpty();  
     rhs.requireWrite();  
 #pragma omp parallel  
     {  
         if (m_faceOffset[0] > -1) {  
             for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring  
 #pragma omp for nowait  
                 for (index_t k1=k1_0; k1<m_NE1; k1+=2) {  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     const index_t e = k1;  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  
                         if (d.actsExpanded()) {  
                             const double d_0 = d_p[0];  
                             const double d_1 = d_p[1];  
                             const double tmp0_0 = d_0 + d_1;  
                             const double tmp1_1 = d_1*w1;  
                             const double tmp4_1 = d_0*w1;  
                             const double tmp0_1 = d_0*w0;  
                             const double tmp3_1 = d_1*w0;  
                             const double tmp2_1 = tmp0_0*w2;  
                             EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp2_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp2_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;  
                         } else { /* constant data */  
                             const double d_0 = d_p[0];  
                             const double tmp1_1 = d_0*w4;  
                             const double tmp0_1 = d_0*w3;  
                             EM_S[INDEX2(0,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp1_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp0_1;  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);  
                         if (y.actsExpanded()) {  
                             const double y_0 = y_p[0];  
                             const double y_1 = y_p[1];  
                             const double tmp3_1 = w5*y_1;  
                             const double tmp2_1 = w6*y_0;  
                             const double tmp0_1 = w6*y_1;  
                             const double tmp1_1 = w5*y_0;  
                             EM_F[0]+=tmp0_1 + tmp1_1;  
                             EM_F[2]+=tmp2_1 + tmp3_1;  
                         } else { /* constant data */  
                             const double y_0 = y_p[0];  
                             const double tmp0_1 = w7*y_0;  
                             EM_F[0]+=tmp0_1;  
                             EM_F[2]+=tmp0_1;  
                         }  
                     }  
                     const index_t firstNode=m_N0*k1;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 }  
             } // end colouring  
         }  
1927    
1928          if (m_faceOffset[1] > -1) {  1234567
1929              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              would be split into two ranks thus:
1930  #pragma omp for nowait  123(4)  (4)567     [4 being a shared element]
                 for (index_t k1=k1_0; k1<m_NE1; k1+=2) {