/[escript]/trunk/ripley/src/Rectangle.cpp

Diff of /trunk/ripley/src/Rectangle.cpp

Parent Directory | Revision Log | View Patch Patch

-revision 4013 by caltinay,
Thu Oct  4 03:13:27 2012 UTC
+revision 4765 by sshaw,
Wed Mar 19 00:17:16 2014 UTC
 Line 1
  /*****************************************************************************
  *
- * Copyright (c) 2003-2012 by University of Queensland
+ * Copyright (c) 2003-2014 by University of Queensland
  * http://www.uq.edu.au
  *
  * Primary Business: Queensland, Australia
 Line 9
  * http://www.opensource.org/licenses/osl-3.0.php
  *
  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
- * Development since 2012 by School of Earth Sciences
+ * Development 2012-2013 by School of Earth Sciences
+ * Development from 2014 by Centre for Geoscience Computing (GeoComp)
  *
  *****************************************************************************/
+ #include <algorithm>
  #include <ripley/Rectangle.h>
- extern "C" {
  #include <paso/SystemMatrix.h>
- }
+ #include <esysUtils/esysFileWriter.h>
+ #include <ripley/DefaultAssembler2D.h>
+ #include <ripley/WaveAssembler2D.h>
+ #include <ripley/LameAssembler2D.h>
+ #include <ripley/domainhelpers.h>
+ #include <boost/scoped_array.hpp>
+ #include "esysUtils/EsysRandom.h"
+ #include "blocktools.h"
  #ifdef USE_NETCDF
  #include <netcdfcpp.h>
-Line 32 
 extern "C" {
+Line 41 
 extern "C" {
  #include <iomanip>
  using namespace std;
+ using esysUtils::FileWriter;
  namespace ripley {
  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
-                      double y1, int d0, int d1) :
+                      double y1, int d0, int d1,
-     RipleyDomain(2),
+                      const std::vector<double>& points,
-     m_x0(x0),
+                      const std::vector<int>& tags,
-     m_y0(y0),
+                      const simap_t& tagnamestonums) :
-     m_l0(x1-x0),
+     RipleyDomain(2)
-     m_l1(y1-y0)
  {
      // ignore subdivision parameters for serial run
      if (m_mpiInfo->size == 1) {
-Line 50 
 Rectangle::Rectangle(int n0, int n1, dou
+Line 59 
 Rectangle::Rectangle(int n0, int n1, dou
      }
      bool warn=false;
-     // if number of subdivisions is non-positive, try to subdivide by the same
+     std::vector<int> factors;
-     // ratio as the number of elements
+     int ranks = m_mpiInfo->size;
-     if (d0<=0 && d1<=0) {
+     int epr[2] = {n0,n1};
-         warn=true;
+     int d[2] = {d0,d1};
-         d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
+     if (d0<=0 || d1<=0) {
-         d1=m_mpiInfo->size/d0;
+         for (int i = 0; i < 2; i++) {
-         if (d0*d1 != m_mpiInfo->size) {
+             if (d[i] < 1) {
-             // ratios not the same so subdivide side with more elements only
+                 d[i] = 1;
-             if (n0>n1) {
+                 continue;
-                 d0=0;
-                 d1=1;
-             } else {
-                 d0=1;
-                 d1=0;
              }
-         }
+             epr[i] = -1; // can no longer be max
-     }
+             //remove
-     if (d0<=0) {
+             if (ranks % d[i] != 0) {
-         // d1 is preset, determine d0 - throw further down if result is no good
+                 throw RipleyException("Invalid number of spatial subdivisions");
-         d0=m_mpiInfo->size/d1;
+             }
-     } else if (d1<=0) {
+             ranks /= d[i];
-         // d0 is preset, determine d1 - throw further down if result is no good
+         }
-         d1=m_mpiInfo->size/d0;
+         factorise(factors, ranks);
+         if (factors.size() != 0) {
+             warn = true;
+         }
+     }
+     while (factors.size() > 0) {
+         int i = epr[0] > epr[1] ? 0 : 1;
+         int f = factors.back();
+         factors.pop_back();
+         d[i] *= f;
+         epr[i] /= f;
      }
+     d0 = d[0]; d1 = d[1];
-     m_NX=d0;
-     m_NY=d1;
      // ensure number of subdivisions is valid and nodes can be distributed
      // among number of ranks
-     if (m_NX*m_NY != m_mpiInfo->size)
+     if (d0*d1 != m_mpiInfo->size)
          throw RipleyException("Invalid number of spatial subdivisions");
      if (warn) {
-Line 88 
 Rectangle::Rectangle(int n0, int n1, dou
+Line 101 
 Rectangle::Rectangle(int n0, int n1, dou
              << d1 << "). This may not be optimal!" << endl;
      }
-     if ((n0+1)%m_NX > 0) {
+     double l0 = x1-x0;
-         double Dx=m_l0/n0;
+     double l1 = y1-y0;
+     m_dx[0] = l0/n0;
+     m_dx[1] = l1/n1;
+     if ((n0+1)%d0 > 0) {
          n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
-         m_l0=Dx*n0;
+         l0=m_dx[0]*n0;
          cout << "Warning: Adjusted number of elements and length. N0="
-             << n0 << ", l0=" << m_l0 << endl;
+             << n0 << ", l0=" << l0 << endl;
      }
-     if ((n1+1)%m_NY > 0) {
+     if ((n1+1)%d1 > 0) {
-         double Dy=m_l1/n1;
          n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
-         m_l1=Dy*n1;
+         l1=m_dx[1]*n1;
          cout << "Warning: Adjusted number of elements and length. N1="
-             << n1 << ", l1=" << m_l1 << endl;
+             << n1 << ", l1=" << l1 << endl;
      }
-     m_gNE0=n0;
+     if ((d0 > 1 && (n0+1)/d0<2) || (d1 > 1 && (n1+1)/d1<2))
-     m_gNE1=n1;
-     if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
          throw RipleyException("Too few elements for the number of ranks");
+     m_gNE[0] = n0;
+     m_gNE[1] = n1;
+     m_origin[0] = x0;
+     m_origin[1] = y0;
+     m_length[0] = l0;
+     m_length[1] = l1;
+     m_NX[0] = d0;
+     m_NX[1] = d1;
      // local number of elements (with and without overlap)
-     m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
+     m_NE[0] = m_ownNE[0] = (d0>1 ? (n0+1)/d0 : n0);
-     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)
-         m_NE0++;
+         m_NE[0]++;
-     else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)
+     else if (d0>1 && m_mpiInfo->rank%d0==d0-1)
-         m_ownNE0--;
+         m_ownNE[0]--;
-     m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
+     m_NE[1] = m_ownNE[1] = (d1>1 ? (n1+1)/d1 : n1);
-     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)
-         m_NE1++;
+         m_NE[1]++;
-     else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)
+     else if (d1>1 && m_mpiInfo->rank/d0==d1-1)
-         m_ownNE1--;
+         m_ownNE[1]--;
      // local number of nodes
-     m_N0 = m_NE0+1;
+     m_NN[0] = m_NE[0]+1;
-     m_N1 = m_NE1+1;
+     m_NN[1] = m_NE[1]+1;
      // bottom-left node is at (offset0,offset1) in global mesh
-     m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);
+     m_offset[0] = (n0+1)/d0*(m_mpiInfo->rank%d0);
-     if (m_offset0 > 0)
+     if (m_offset[0] > 0)
-         m_offset0--;
+         m_offset[0]--;
-     m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);
+     m_offset[1] = (n1+1)/d1*(m_mpiInfo->rank/d0);
-     if (m_offset1 > 0)
+     if (m_offset[1] > 0)
-         m_offset1--;
+         m_offset[1]--;
      populateSampleIds();
      createPattern();
+     assembler = new DefaultAssembler2D(this, m_dx, m_NX, m_NE, m_NN);
+     for (map<string, int>::const_iterator i = tagnamestonums.begin();
+             i != tagnamestonums.end(); i++) {
+         setTagMap(i->first, i->second);
+     }
+     addPoints(tags.size(), &points[0], &tags[0]);
  }
  Rectangle::~Rectangle()
  {
      Paso_SystemMatrixPattern_free(m_pattern);
      Paso_Connector_free(m_connector);
+     delete assembler;
  }
  string Rectangle::getDescription() const
-Line 154 
 bool Rectangle::operator==(const Abstrac
+Line 183 
 bool Rectangle::operator==(const Abstrac
      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
      if (o) {
          return (RipleyDomain::operator==(other) &&
-                 m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1
+                 m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1]
-                 && m_x0==o->m_x0 && m_y0==o->m_y0
+                 && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1]
-                 && m_l0==o->m_l0 && m_l1==o->m_l1
+                 && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1]
-                 && m_NX==o->m_NX && m_NY==o->m_NY);
+                 && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1]);
      }
      return false;
  }
  void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
-             const vector<int>& first, const vector<int>& numValues) const
+             const ReaderParameters& params) const
  {
  #ifdef USE_NETCDF
      // check destination function space
      int myN0, myN1;
      if (out.getFunctionSpace().getTypeCode() == Nodes) {
-         myN0 = m_N0;
+         myN0 = m_NN[0];
-         myN1 = m_N1;
+         myN1 = m_NN[1];
      } else if (out.getFunctionSpace().getTypeCode() == Elements ||
                  out.getFunctionSpace().getTypeCode() == ReducedElements) {
-         myN0 = m_NE0;
+         myN0 = m_NE[0];
-         myN1 = m_NE1;
+         myN1 = m_NE[1];
      } else
          throw RipleyException("readNcGrid(): invalid function space for output data object");
-     if (first.size() != 2)
+     if (params.first.size() != 2)
          throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
-     if (numValues.size() != 2)
+     if (params.numValues.size() != 2)
          throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
+     if (params.multiplier.size() != 2)
+         throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");
+     for (size_t i=0; i<params.multiplier.size(); i++)
+         if (params.multiplier[i]<1)
+             throw RipleyException("readNcGrid(): all multipliers must be positive");
+     if (params.reverse.size() != 2)
+         throw RipleyException("readNcGrid(): argument 'reverse' must have 2 entries");
      // check file existence and size
      NcFile f(filename.c_str(), NcFile::ReadOnly);
      if (!f.is_valid())
-Line 200 
 void Rectangle::readNcGrid(escript::Data
+Line 237 
 void Rectangle::readNcGrid(escript::Data
          throw RipleyException("readNcGrid(): only scalar data supported");
      const int dims = var->num_dims();
-     const long *edges = var->edges();
+     boost::scoped_array<long> edges(var->edges());
      // is this a slice of the data object (dims!=2)?
      // note the expected ordering of edges (as in numpy: y,x)
-     if ( (dims==2 && (numValues[1] > edges[0] || numValues[0] > edges[1]))
+     if ( (dims==2 && (params.numValues[1] > edges[0] || params.numValues[0] > edges[1]))
-             || (dims==1 && numValues[1]>1) ) {
+             || (dims==1 && params.numValues[1]>1) ) {
          throw RipleyException("readNcGrid(): not enough data in file");
      }
      // check if this rank contributes anything
-     if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
+     if (params.first[0] >= m_offset[0]+myN0 ||
-             first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1)
+             params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
+             params.first[1] >= m_offset[1]+myN1 ||
+             params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1])
          return;
      // now determine how much this rank has to write
      // first coordinates in data object to write to
-     const int first0 = max(0, first[0]-m_offset0);
+     const int first0 = max(0, params.first[0]-m_offset[0]);
-     const int first1 = max(0, first[1]-m_offset1);
+     const int first1 = max(0, params.first[1]-m_offset[1]);
-     // indices to first value in file
+     // indices to first value in file (not accounting for reverse yet)
-     const int idx0 = max(0, m_offset0-first[0]);
+     int idx0 = max(0, m_offset[0]-params.first[0]);
-     const int idx1 = max(0, m_offset1-first[1]);
+     int idx1 = max(0, m_offset[1]-params.first[1]);
-     // number of values to write
+     // number of values to read
-     const int num0 = min(numValues[0]-idx0, myN0-first0);
+     const int num0 = min(params.numValues[0]-idx0, myN0-first0);
-     const int num1 = min(numValues[1]-idx1, myN1-first1);
+     const int num1 = min(params.numValues[1]-idx1, myN1-first1);
+     // make sure we read the right block if going backwards through file
+     if (params.reverse[0])
+         idx0 = edges[dims-1]-num0-idx0;
+     if (dims>1 && params.reverse[1])
+         idx1 = edges[dims-2]-num1-idx1;
      vector<double> values(num0*num1);
      if (dims==2) {
-Line 238 
 void Rectangle::readNcGrid(escript::Data
+Line 283 
 void Rectangle::readNcGrid(escript::Data
      const int dpp = out.getNumDataPointsPerSample();
      out.requireWrite();
+     // helpers for reversing
+     const int x0 = (params.reverse[0] ? num0-1 : 0);
+     const int x_mult = (params.reverse[0] ? -1 : 1);
+     const int y0 = (params.reverse[1] ? num1-1 : 0);
+     const int y_mult = (params.reverse[1] ? -1 : 1);
      for (index_t y=0; y<num1; y++) {
  #pragma omp parallel for
          for (index_t x=0; x<num0; x++) {
-             const int dataIndex = (first1+y)*myN0+first0+x;
+             const int baseIndex = first0+x*params.multiplier[0]
-             const int srcIndex=y*num0+x;
+                                   +(first1+y*params.multiplier[1])*myN0;
-             double* dest = out.getSampleDataRW(dataIndex);
+             const int srcIndex = (y0+y_mult*y)*num0+(x0+x_mult*x);
-             for (index_t q=0; q<dpp; q++) {
+             if (!isnan(values[srcIndex])) {
-                 *dest++ = values[srcIndex];
+                 for (index_t m1=0; m1<params.multiplier[1]; m1++) {
+                     for (index_t m0=0; m0<params.multiplier[0]; m0++) {
+                         const int dataIndex = baseIndex+m0+m1*myN0;
+                         double* dest = out.getSampleDataRW(dataIndex);
+                         for (index_t q=0; q<dpp; q++) {
+                             *dest++ = values[srcIndex];
+                         }
+                     }
+                 }
              }
          }
      }
-Line 255 
 void Rectangle::readNcGrid(escript::Data
+Line 314 
 void Rectangle::readNcGrid(escript::Data
  }
  void Rectangle::readBinaryGrid(escript::Data& out, string filename,
-             const vector<int>& first, const vector<int>& numValues) const
+                                const ReaderParameters& params) const
+ {
+     // the mapping is not universally correct but should work on our
+     // supported platforms
+     switch (params.dataType) {
+         case DATATYPE_INT32:
+             readBinaryGridImpl<int>(out, filename, params);
+             break;
+         case DATATYPE_FLOAT32:
+             readBinaryGridImpl<float>(out, filename, params);
+             break;
+         case DATATYPE_FLOAT64:
+             readBinaryGridImpl<double>(out, filename, params);
+             break;
+         default:
+             throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
+     }
+ }
+ void Rectangle::readBinaryGridFromZipped(escript::Data& out, string filename,
+                                const ReaderParameters& params) const
+ {
+     // the mapping is not universally correct but should work on our
+     // supported platforms
+     switch (params.dataType) {
+         case DATATYPE_INT32:
+             readBinaryGridZippedImpl<int>(out, filename, params);
+             break;
+         case DATATYPE_FLOAT32:
+             readBinaryGridZippedImpl<float>(out, filename, params);
+             break;
+         case DATATYPE_FLOAT64:
+             readBinaryGridZippedImpl<double>(out, filename, params);
+             break;
+         default:
+             throw RipleyException("readBinaryGridFromZipped(): invalid or unsupported datatype");
+     }
+ }
+ template<typename ValueType>
+ void Rectangle::readBinaryGridImpl(escript::Data& out, const string& filename,
+                                    const ReaderParameters& params) const
  {
      // check destination function space
      int myN0, myN1;
      if (out.getFunctionSpace().getTypeCode() == Nodes) {
-         myN0 = m_N0;
+         myN0 = m_NN[0];
-         myN1 = m_N1;
+         myN1 = m_NN[1];
      } else if (out.getFunctionSpace().getTypeCode() == Elements ||
                  out.getFunctionSpace().getTypeCode() == ReducedElements) {
-         myN0 = m_NE0;
+         myN0 = m_NE[0];
-         myN1 = m_NE1;
+         myN1 = m_NE[1];
      } else
          throw RipleyException("readBinaryGrid(): invalid function space for output data object");
-Line 277 
 void Rectangle::readBinaryGrid(escript::
+Line 377 
 void Rectangle::readBinaryGrid(escript::
      f.seekg(0, ios::end);
      const int numComp = out.getDataPointSize();
      const int filesize = f.tellg();
-     const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(float);
+     const int reqsize = params.numValues[0]*params.numValues[1]*numComp*sizeof(ValueType);
      if (filesize < reqsize) {
          f.close();
          throw RipleyException("readBinaryGrid(): not enough data in file");
      }
      // check if this rank contributes anything
-     if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
+     if (params.first[0] >= m_offset[0]+myN0 ||
-             first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1) {
+             params.first[0]+params.numValues[0] <= m_offset[0] ||
+             params.first[1] >= m_offset[1]+myN1 ||
+             params.first[1]+params.numValues[1] <= m_offset[1]) {
          f.close();
          return;
      }
-Line 293 
 void Rectangle::readBinaryGrid(escript::
+Line 395 
 void Rectangle::readBinaryGrid(escript::
      // now determine how much this rank has to write
      // first coordinates in data object to write to
-     const int first0 = max(0, first[0]-m_offset0);
+     const int first0 = max(0, params.first[0]-m_offset[0]);
-     const int first1 = max(0, first[1]-m_offset1);
+     const int first1 = max(0, params.first[1]-m_offset[1]);
      // indices to first value in file
-     const int idx0 = max(0, m_offset0-first[0]);
+     const int idx0 = max(0, m_offset[0]-params.first[0]);
-     const int idx1 = max(0, m_offset1-first[1]);
+     const int idx1 = max(0, m_offset[1]-params.first[1]);
-     // number of values to write
+     // number of values to read
-     const int num0 = min(numValues[0]-idx0, myN0-first0);
+     const int num0 = min(params.numValues[0]-idx0, myN0-first0);
-     const int num1 = min(numValues[1]-idx1, myN1-first1);
+     const int num1 = min(params.numValues[1]-idx1, myN1-first1);
      out.requireWrite();
-     vector<float> values(num0*numComp);
+     vector<ValueType> values(num0*numComp);
      const int dpp = out.getNumDataPointsPerSample();
-     for (index_t y=0; y<num1; y++) {
+     for (int y=0; y<num1; y++) {
-         const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);
+         const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]);
-         f.seekg(fileofs*sizeof(float));
+         f.seekg(fileofs*sizeof(ValueType));
-         f.read((char*)&values[0], num0*numComp*sizeof(float));
+         f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
-         for (index_t x=0; x<num0; x++) {
+         for (int x=0; x<num0; x++) {
-             double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0);
+             const int baseIndex = first0+x*params.multiplier[0]
-             for (index_t c=0; c<numComp; c++) {
+                                     +(first1+y*params.multiplier[1])*myN0;
-                 for (index_t q=0; q<dpp; q++) {
+             for (int m1=0; m1<params.multiplier[1]; m1++) {
-                     *dest++ = static_cast<double>(values[x*numComp+c]);
+                 for (int m0=0; m0<params.multiplier[0]; m0++) {
+                     const int dataIndex = baseIndex+m0+m1*myN0;
+                     double* dest = out.getSampleDataRW(dataIndex);
+                     for (int c=0; c<numComp; c++) {
+                         ValueType val = values[x*numComp+c];
+                         if (params.byteOrder != BYTEORDER_NATIVE) {
+                             char* cval = reinterpret_cast<char*>(&val);
+                             // this will alter val!!
+                             byte_swap32(cval);
+                         }
+                         if (!std::isnan(val)) {
+                             for (int q=0; q<dpp; q++) {
+                                 *dest++ = static_cast<double>(val);
+                             }
+                         }
+                     }
+                 }
+             }
+         }
+     }
+     f.close();
+ }
+ template<typename ValueType>
+ void Rectangle::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
+                                    const ReaderParameters& params) const
+ {
+     // check destination function space
+     int myN0, myN1;
+     if (out.getFunctionSpace().getTypeCode() == Nodes) {
+         myN0 = m_NN[0];
+         myN1 = m_NN[1];
+     } else if (out.getFunctionSpace().getTypeCode() == Elements ||
+                 out.getFunctionSpace().getTypeCode() == ReducedElements) {
+         myN0 = m_NE[0];
+         myN1 = m_NE[1];
+     } else
+         throw RipleyException("readBinaryGrid(): invalid function space for output data object");
+     // check file existence and size
+     ifstream f(filename.c_str(), ifstream::binary);
+     if (f.fail()) {
+         throw RipleyException("readBinaryGridFromZipped(): cannot open file");
+     }
+     f.seekg(0, ios::end);
+     const int numComp = out.getDataPointSize();
+     int filesize = f.tellg();
+     f.seekg(0, ios::beg);
+     std::vector<char> compressed(filesize);
+     f.read((char*)&compressed[0], filesize);
+     f.close();
+     std::vector<char> decompressed = unzip(compressed);
+     filesize = decompressed.size();
+     const int reqsize = params.numValues[0]*params.numValues[1]*numComp*sizeof(ValueType);
+     if (filesize < reqsize) {
+         throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
+     }
+     // check if this rank contributes anything
+     if (params.first[0] >= m_offset[0]+myN0 ||
+             params.first[0]+params.numValues[0] <= m_offset[0] ||
+             params.first[1] >= m_offset[1]+myN1 ||
+             params.first[1]+params.numValues[1] <= m_offset[1]) {
+         f.close();
+         return;
+     }
+     // now determine how much this rank has to write
+     // first coordinates in data object to write to
+     const int first0 = max(0, params.first[0]-m_offset[0]);
+     const int first1 = max(0, params.first[1]-m_offset[1]);
+     // indices to first value in file
+     const int idx0 = max(0, m_offset[0]-params.first[0]);
+     const int idx1 = max(0, m_offset[1]-params.first[1]);
+     // number of values to read
+     const int num0 = min(params.numValues[0]-idx0, myN0-first0);
+     const int num1 = min(params.numValues[1]-idx1, myN1-first1);
+     out.requireWrite();
+     vector<ValueType> values(num0*numComp);
+     const int dpp = out.getNumDataPointsPerSample();
+     for (int y=0; y<num1; y++) {
+         const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]);
+             memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
+         for (int x=0; x<num0; x++) {
+             const int baseIndex = first0+x*params.multiplier[0]
+                                     +(first1+y*params.multiplier[1])*myN0;
+             for (int m1=0; m1<params.multiplier[1]; m1++) {
+                 for (int m0=0; m0<params.multiplier[0]; m0++) {
+                     const int dataIndex = baseIndex+m0+m1*myN0;
+                     double* dest = out.getSampleDataRW(dataIndex);
+                     for (int c=0; c<numComp; c++) {
+                         ValueType val = values[x*numComp+c];
+                         if (params.byteOrder != BYTEORDER_NATIVE) {
+                             char* cval = reinterpret_cast<char*>(&val);
+                             // this will alter val!!
+                             byte_swap32(cval);
+                         }
+                         if (!std::isnan(val)) {
+                             for (int q=0; q<dpp; q++) {
+                                 *dest++ = static_cast<double>(val);
+                             }
+                         }
+                     }
                  }
              }
          }
-Line 323 
 void Rectangle::readBinaryGrid(escript::
+Line 533 
 void Rectangle::readBinaryGrid(escript::
      f.close();
  }
+ void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,
+                                 int byteOrder, int dataType) const
+ {
+     // the mapping is not universally correct but should work on our
+     // supported platforms
+     switch (dataType) {
+         case DATATYPE_INT32:
+             writeBinaryGridImpl<int>(in, filename, byteOrder);
+             break;
+         case DATATYPE_FLOAT32:
+             writeBinaryGridImpl<float>(in, filename, byteOrder);
+             break;
+         case DATATYPE_FLOAT64:
+             writeBinaryGridImpl<double>(in, filename, byteOrder);
+             break;
+         default:
+             throw RipleyException("writeBinaryGrid(): invalid or unsupported datatype");
+     }
+ }
+ template<typename ValueType>
+ void Rectangle::writeBinaryGridImpl(const escript::Data& in,
+                                     const string& filename, int byteOrder) const
+ {
+     // check function space and determine number of points
+     int myN0, myN1;
+     int totalN0, totalN1;
+     if (in.getFunctionSpace().getTypeCode() == Nodes) {
+         myN0 = m_NN[0];
+         myN1 = m_NN[1];
+         totalN0 = m_gNE[0]+1;
+         totalN1 = m_gNE[1]+1;
+     } else if (in.getFunctionSpace().getTypeCode() == Elements ||
+                 in.getFunctionSpace().getTypeCode() == ReducedElements) {
+         myN0 = m_NE[0];
+         myN1 = m_NE[1];
+         totalN0 = m_gNE[0];
+         totalN1 = m_gNE[1];
+     } else
+         throw RipleyException("writeBinaryGrid(): invalid function space of data object");
+     const int numComp = in.getDataPointSize();
+     const int dpp = in.getNumDataPointsPerSample();
+     if (numComp > 1 || dpp > 1)
+         throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
+     const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;
+     // from here on we know that each sample consists of one value
+     FileWriter fw;
+     fw.openFile(filename, fileSize);
+     MPIBarrier();
+     for (index_t y=0; y<myN1; y++) {
+         const int fileofs = (m_offset[0]+(m_offset[1]+y)*totalN0)*sizeof(ValueType);
+         ostringstream oss;
+         for (index_t x=0; x<myN0; x++) {
+             const double* sample = in.getSampleDataRO(y*myN0+x);
+             ValueType fvalue = static_cast<ValueType>(*sample);
+             if (byteOrder == BYTEORDER_NATIVE) {
+                 oss.write((char*)&fvalue, sizeof(fvalue));
+             } else {
+                 char* value = reinterpret_cast<char*>(&fvalue);
+                 oss.write(byte_swap32(value), sizeof(fvalue));
+             }
+         }
+         fw.writeAt(oss, fileofs);
+     }
+     fw.close();
+ }
  void Rectangle::dump(const string& fileName) const
  {
  #if USE_SILO
-Line 385 
 void Rectangle::dump(const string& fileN
+Line 668 
 void Rectangle::dump(const string& fileN
      }
      */
-     boost::scoped_ptr<double> x(new double[m_N0]);
+     boost::scoped_ptr<double> x(new double[m_NN[0]]);
-     boost::scoped_ptr<double> y(new double[m_N1]);
+     boost::scoped_ptr<double> y(new double[m_NN[1]]);
      double* coords[2] = { x.get(), y.get() };
-     pair<double,double> xdx = getFirstCoordAndSpacing(0);
-     pair<double,double> ydy = getFirstCoordAndSpacing(1);
  #pragma omp parallel
      {
  #pragma omp for nowait
-         for (dim_t i0 = 0; i0 < m_N0; i0++) {
+         for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
-             coords[0][i0]=xdx.first+i0*xdx.second;
+             coords[0][i0]=getLocalCoordinate(i0, 0);
          }
  #pragma omp for nowait
-         for (dim_t i1 = 0; i1 < m_N1; i1++) {
+         for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
-             coords[1][i1]=ydy.first+i1*ydy.second;
+             coords[1][i1]=getLocalCoordinate(i1, 1);
          }
      }
-     IndexVector dims = getNumNodesPerDim();
+     int* dims = const_cast<int*>(getNumNodesPerDim());
      // write mesh
-     DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,
+     DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,
              DB_COLLINEAR, NULL);
      // write node ids
-     DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 2,
+     DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 2,
              NULL, 0, DB_INT, DB_NODECENT, NULL);
      // write element ids
-     dims = getNumElementsPerDim();
+     dims = const_cast<int*>(getNumElementsPerDim());
      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
-             &dims[0], 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);
+             dims, 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);
      // rank 0 writes multimesh and multivar
      if (m_mpiInfo->rank == 0) {
-Line 482 
 const int* Rectangle::borrowSampleRefere
+Line 763 
 const int* Rectangle::borrowSampleRefere
          case FaceElements:
          case ReducedFaceElements:
              return &m_faceId[0];
+         case Points:
+             return &m_diracPointNodeIDs[0];
          default:
              break;
      }
-Line 506 
 bool Rectangle::ownSample(int fsType, in
+Line 789 
 bool Rectangle::ownSample(int fsType, in
          case Elements:
          case ReducedElements:
              // check ownership of element's bottom left node
-             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());
          case FaceElements:
          case ReducedFaceElements:
              {
                  // determine which face the sample belongs to before
                  // checking ownership of corresponding element's first node
-                 const IndexVector faces = getNumFacesPerBoundary();
                  dim_t n=0;
-                 for (size_t i=0; i<faces.size(); i++) {
+                 for (size_t i=0; i<4; i++) {
-                     n+=faces[i];
+                     n+=m_faceCount[i];
                      if (id<n) {
                          index_t k;
                          if (i==1)
-                             k=m_N0-2;
+                             k=m_NN[0]-2;
                          else if (i==3)
-                             k=m_N0*(m_N1-2);
+                             k=m_NN[0]*(m_NN[1]-2);
                          else
                              k=0;
                          // determine whether to move right or up
-                         const index_t delta=(i/2==0 ? m_N0 : 1);
+                         const index_t delta=(i/2==0 ? m_NN[0] : 1);
-                         return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());
+                         return (m_dofMap[k+(id-n+m_faceCount[i])*delta] < getNumDOF());
                      }
                  }
                  return false;
-Line 548 
 void Rectangle::setToNormal(escript::Dat
+Line 830 
 void Rectangle::setToNormal(escript::Dat
          {
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
+                 for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
                      // set vector at two quadrature points
                      *o++ = -1.;
-Line 560 
 void Rectangle::setToNormal(escript::Dat
+Line 842 
 void Rectangle::setToNormal(escript::Dat
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
+                 for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
                      // set vector at two quadrature points
                      *o++ = 1.;
-Line 572 
 void Rectangle::setToNormal(escript::Dat
+Line 854 
 void Rectangle::setToNormal(escript::Dat
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
+                 for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
                      // set vector at two quadrature points
                      *o++ = 0.;
-Line 584 
 void Rectangle::setToNormal(escript::Dat
+Line 866 
 void Rectangle::setToNormal(escript::Dat
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
+                 for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
                      // set vector at two quadrature points
                      *o++ = 0.;
-Line 600 
 void Rectangle::setToNormal(escript::Dat
+Line 882 
 void Rectangle::setToNormal(escript::Dat
          {
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
+                 for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
                      *o++ = -1.;
                      *o = 0.;
-Line 609 
 void Rectangle::setToNormal(escript::Dat
+Line 891 
 void Rectangle::setToNormal(escript::Dat
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
+                 for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
                      *o++ = 1.;
                      *o = 0.;
-Line 618 
 void Rectangle::setToNormal(escript::Dat
+Line 900 
 void Rectangle::setToNormal(escript::Dat
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
+                 for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
                      *o++ = 0.;
                      *o = -1.;
-Line 627 
 void Rectangle::setToNormal(escript::Dat
+Line 909 
 void Rectangle::setToNormal(escript::Dat
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
+                 for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
                      *o++ = 0.;
                      *o = 1.;
-Line 649 
 void Rectangle::setToSize(escript::Data&
+Line 931 
 void Rectangle::setToSize(escript::Data&
              || out.getFunctionSpace().getTypeCode() == ReducedElements) {
          out.requireWrite();
          const dim_t numQuad=out.getNumDataPointsPerSample();
-         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=sqrt(hSize*hSize+vSize*vSize);
  #pragma omp parallel for
          for (index_t k = 0; k < getNumElements(); ++k) {
              double* o = out.getSampleDataRW(k);
-Line 661 
 void Rectangle::setToSize(escript::Data&
+Line 941 
 void Rectangle::setToSize(escript::Data&
              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
          out.requireWrite();
          const dim_t numQuad=out.getNumDataPointsPerSample();
-         const double hSize=getFirstCoordAndSpacing(0).second;
-         const double vSize=getFirstCoordAndSpacing(1).second;
  #pragma omp parallel
          {
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
+                 for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
-                     fill(o, o+numQuad, vSize);
+                     fill(o, o+numQuad, m_dx[1]);
                  }
              }
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
+                 for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
-                     fill(o, o+numQuad, vSize);
+                     fill(o, o+numQuad, m_dx[1]);
                  }
              }
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
+                 for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
-                     fill(o, o+numQuad, hSize);
+                     fill(o, o+numQuad, m_dx[0]);
                  }
              }
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
+                 for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
-                     fill(o, o+numQuad, hSize);
+                     fill(o, o+numQuad, m_dx[0]);
                  }
              }
          } // end of parallel section
-Line 706 
 void Rectangle::setToSize(escript::Data&
+Line 984 
 void Rectangle::setToSize(escript::Data&
      }
  }
- 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;
- }
  void Rectangle::Print_Mesh_Info(const bool full) const
  {
      RipleyDomain::Print_Mesh_Info(full);
-Line 723 
 void Rectangle::Print_Mesh_Info(const bo
+Line 991 
 void Rectangle::Print_Mesh_Info(const bo
          cout << "     Id  Coordinates" << endl;
          cout.precision(15);
          cout.setf(ios::scientific, ios::floatfield);
-         pair<double,double> xdx = getFirstCoordAndSpacing(0);
-         pair<double,double> ydy = getFirstCoordAndSpacing(1);
          for (index_t i=0; i < getNumNodes(); i++) {
              cout << "  " << setw(5) << m_nodeId[i]
-                 << "  " << xdx.first+(i%m_N0)*xdx.second
+                 << "  " << getLocalCoordinate(i%m_NN[0], 0)
-                 << "  " << ydy.first+(i/m_N0)*ydy.second << endl;
+                 << "  " << getLocalCoordinate(i/m_NN[0], 1) << endl;
          }
      }
  }
- 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;
- }
  //protected
  void Rectangle::assembleCoordinates(escript::Data& arg) const
-Line 811 
 void Rectangle::assembleCoordinates(escr
+Line 1010 
 void Rectangle::assembleCoordinates(escr
      if (!numSamplesEqual(&x, 1, getNumNodes()))
          throw RipleyException("setToX: Illegal number of samples in Data object");
-     pair<double,double> xdx = getFirstCoordAndSpacing(0);
-     pair<double,double> ydy = getFirstCoordAndSpacing(1);
      arg.requireWrite();
  #pragma omp parallel for
-     for (dim_t i1 = 0; i1 < m_N1; i1++) {
+     for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
-         for (dim_t i0 = 0; i0 < m_N0; i0++) {
+         for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
-             double* point = arg.getSampleDataRW(i0+m_N0*i1);
+             double* point = arg.getSampleDataRW(i0+m_NN[0]*i1);
-             point[0] = xdx.first+i0*xdx.second;
+             point[0] = getLocalCoordinate(i0, 0);
-             point[1] = ydy.first+i1*ydy.second;
+             point[1] = getLocalCoordinate(i1, 1);
          }
      }
  }
  //protected
- void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
+ void Rectangle::assembleGradient(escript::Data& out, const escript::Data& in) const
  {
      const dim_t numComp = in.getDataPointSize();
-     const double h0 = m_l0/m_gNE0;
+     const double cx0 = .21132486540518711775/m_dx[0];
-     const double h1 = m_l1/m_gNE1;
+     const double cx1 = .78867513459481288225/m_dx[0];
-     const double cx0 = -1./h0;
+     const double cx2 = 1./m_dx[0];
-     const double cx1 = -.78867513459481288225/h0;
+     const double cy0 = .21132486540518711775/m_dx[1];
-     const double cx2 = -.5/h0;
+     const double cy1 = .78867513459481288225/m_dx[1];
-     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;
      if (out.getFunctionSpace().getTypeCode() == Elements) {
          out.requireWrite();
-Line 856 
 void Rectangle::assembleGradient(escript
+Line 1041 
 void Rectangle::assembleGradient(escript
              vector<double> f_10(numComp);
              vector<double> f_11(numComp);
  #pragma omp for
-             for (index_t k1=0; k1 < m_NE1; ++k1) {
+             for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
-                     double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
+                     double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
-                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
+                         o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
-                         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;
-                         o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
+                         o[INDEX3(i,0,2,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
-                         o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
+                         o[INDEX3(i,1,2,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
-                         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,3,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
-                         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,3,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
                      } // end of component loop i
                  } // end of k0 loop
              } // end of k1 loop
-Line 885 
 void Rectangle::assembleGradient(escript
+Line 1070 
 void Rectangle::assembleGradient(escript
              vector<double> f_10(numComp);
              vector<double> f_11(numComp);
  #pragma omp for
-             for (index_t k1=0; k1 < m_NE1; ++k1) {
+             for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
-                     double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
+                     double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
                      } // end of component loop i
                  } // end of k0 loop
              } // end of k1 loop
-Line 909 
 void Rectangle::assembleGradient(escript
+Line 1094 
 void Rectangle::assembleGradient(escript
              vector<double> f_11(numComp);
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                 for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
-                         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;
-                         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;
                      } // end of component loop i
                  } // end of k1 loop
              } // end of face 0
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                 for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
-                         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;
-                         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;
                      } // end of component loop i
                  } // end of k1 loop
              } // end of face 1
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
-                         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;
-                         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;
                      } // end of component loop i
                  } // end of k0 loop
              } // end of face 2
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
-                         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;
-                         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;
                      } // end of component loop i
                  } // end of k0 loop
              } // end of face 3
-Line 983 
 void Rectangle::assembleGradient(escript
+Line 1168 
 void Rectangle::assembleGradient(escript
              vector<double> f_11(numComp);
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                 for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
                      } // end of component loop i
                  } // end of k1 loop
              } // end of face 0
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                 for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
                      } // end of component loop i
                  } // end of k1 loop
              } // end of face 1
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
                      } // end of component loop i
                  } // end of k0 loop
              } // end of face 2
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
                      for (index_t i=0; i < numComp; ++i) {
-                         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;
-                         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;
                      } // end of component loop i
                  } // end of k0 loop
              } // end of face 3
-Line 1042 
 void Rectangle::assembleGradient(escript
+Line 1227 
 void Rectangle::assembleGradient(escript
  }
  //protected
- void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
+ void Rectangle::assembleIntegrate(vector<double>& integrals,
+                                   const escript::Data& arg) const
  {
      const dim_t numComp = arg.getDataPointSize();
-     const double h0 = m_l0/m_gNE0;
+     const index_t left = (m_offset[0]==0 ? 0 : 1);
-     const double h1 = m_l1/m_gNE1;
+     const index_t bottom = (m_offset[1]==0 ? 0 : 1);
-     const index_t left = (m_offset0==0 ? 0 : 1);
-     const index_t bottom = (m_offset1==0 ? 0 : 1);
      const int fs=arg.getFunctionSpace().getTypeCode();
      if (fs == Elements && arg.actsExpanded()) {
  #pragma omp parallel
          {
              vector<double> int_local(numComp, 0);
-             const double w = h0*h1/4.;
+             const double w = m_dx[0]*m_dx[1]/4.;
  #pragma omp for nowait
-             for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+             for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
-                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                 for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
-                     const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
+                     const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
                      for (index_t i=0; i < numComp; ++i) {
                          const double f0 = f[INDEX2(i,0,numComp)];
                          const double f1 = f[INDEX2(i,1,numComp)];
-Line 1074 
 void Rectangle::assembleIntegrate(vector
+Line 1258 
 void Rectangle::assembleIntegrate(vector
          } // end of parallel section
      } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
-         const double w = h0*h1;
+         const double w = m_dx[0]*m_dx[1];
  #pragma omp parallel
          {
              vector<double> int_local(numComp, 0);
  #pragma omp for nowait
-             for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+             for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
-                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                 for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
-                     const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
+                     const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
                      for (index_t i=0; i < numComp; ++i) {
                          int_local[i]+=f[i]*w;
                      }
-Line 1096 
 void Rectangle::assembleIntegrate(vector
+Line 1280 
 void Rectangle::assembleIntegrate(vector
  #pragma omp parallel
          {
              vector<double> int_local(numComp, 0);
-             const double w0 = h0/2.;
+             const double w0 = m_dx[0]/2.;
-             const double w1 = h1/2.;
+             const double w1 = m_dx[1]/2.;
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                 for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
                      for (index_t i=0; i < numComp; ++i) {
                          const double f0 = f[INDEX2(i,0,numComp)];
-Line 1112 
 void Rectangle::assembleIntegrate(vector
+Line 1296 
 void Rectangle::assembleIntegrate(vector
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                 for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
                      for (index_t i=0; i < numComp; ++i) {
                          const double f0 = f[INDEX2(i,0,numComp)];
-Line 1124 
 void Rectangle::assembleIntegrate(vector
+Line 1308 
 void Rectangle::assembleIntegrate(vector
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                 for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
                      for (index_t i=0; i < numComp; ++i) {
                          const double f0 = f[INDEX2(i,0,numComp)];
-Line 1136 
 void Rectangle::assembleIntegrate(vector
+Line 1320 
 void Rectangle::assembleIntegrate(vector
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                 for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
                      for (index_t i=0; i < numComp; ++i) {
                          const double f0 = f[INDEX2(i,0,numComp)];
-Line 1156 
 void Rectangle::assembleIntegrate(vector
+Line 1340 
 void Rectangle::assembleIntegrate(vector
              vector<double> int_local(numComp, 0);
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                 for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
                      for (index_t i=0; i < numComp; ++i) {
-                         int_local[i]+=f[i]*h1;
+                         int_local[i]+=f[i]*m_dx[1];
                      }
                  }
              }
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
-                 for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                 for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
                      for (index_t i=0; i < numComp; ++i) {
-                         int_local[i]+=f[i]*h1;
+                         int_local[i]+=f[i]*m_dx[1];
                      }
                  }
              }
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                 for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
                      for (index_t i=0; i < numComp; ++i) {
-                         int_local[i]+=f[i]*h0;
+                         int_local[i]+=f[i]*m_dx[0];
                      }
                  }
              }
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
-                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                 for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
                      for (index_t i=0; i < numComp; ++i) {
-                         int_local[i]+=f[i]*h0;
+                         int_local[i]+=f[i]*m_dx[0];
                      }
                  }
              }
-Line 1204 
 void Rectangle::assembleIntegrate(vector
+Line 1388 
 void Rectangle::assembleIntegrate(vector
  //protected
  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
  {
-     const dim_t nDOF0 = (m_gNE0+1)/m_NX;
+     const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
-     const dim_t nDOF1 = (m_gNE1+1)/m_NY;
+     const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
      const int x=node%nDOF0;
      const int y=node/nDOF0;
      dim_t num=0;
-Line 1230 
 dim_t Rectangle::insertNeighbourNodes(In
+Line 1414 
 dim_t Rectangle::insertNeighbourNodes(In
  }
  //protected
- void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const
+ void Rectangle::nodesToDOF(escript::Data& out, const escript::Data& in) const
  {
      const dim_t numComp = in.getDataPointSize();
      out.requireWrite();
-     const index_t left = (m_offset0==0 ? 0 : 1);
+     const index_t left = (m_offset[0]==0 ? 0 : 1);
-     const index_t bottom = (m_offset1==0 ? 0 : 1);
+     const index_t bottom = (m_offset[1]==0 ? 0 : 1);
-     const dim_t nDOF0 = (m_gNE0+1)/m_NX;
+     const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
-     const dim_t nDOF1 = (m_gNE1+1)/m_NY;
+     const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
  #pragma omp parallel for
      for (index_t i=0; i<nDOF1; i++) {
          for (index_t j=0; j<nDOF0; j++) {
-             const index_t n=j+left+(i+bottom)*m_N0;
+             const index_t n=j+left+(i+bottom)*m_NN[0];
              const double* src=in.getSampleDataRO(n);
              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
          }
-Line 1250 
 void Rectangle::nodesToDOF(escript::Data
+Line 1434 
 void Rectangle::nodesToDOF(escript::Data
  }
  //protected
- void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const
+ void Rectangle::dofToNodes(escript::Data& out, const escript::Data& in) const
  {
      const dim_t numComp = in.getDataPointSize();
      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
-     in.requireWrite();
+     // expand data object if necessary to be able to grab the whole data
-     Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));
+     const_cast<escript::Data*>(&in)->expand();
+     Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));
      const dim_t numDOF = getNumDOF();
      out.requireWrite();
-Line 1274 
 void Rectangle::dofToNodes(escript::Data
+Line 1459 
 void Rectangle::dofToNodes(escript::Data
  //private
  void Rectangle::populateSampleIds()
  {
-     // identifiers are ordered from left to right, bottom to top globablly.
+     // degrees of freedom are numbered from left to right, bottom to top in
+     // each rank, continuing on the next rank (ranks also go left-right,
+     // bottom-top).
+     // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
+     // helps when writing out data rank after rank.
      // build node distribution vector first.
-     // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
+     // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
+     // constant for all ranks in this implementation
      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
      const dim_t numDOF=getNumDOF();
      for (dim_t k=1; k<m_mpiInfo->size; k++) {
-Line 1287 
 void Rectangle::populateSampleIds()
+Line 1477 
 void Rectangle::populateSampleIds()
      m_nodeId.resize(getNumNodes());
      m_dofId.resize(numDOF);
      m_elementId.resize(getNumElements());
+     // populate face element counts
+     //left
+     if (m_offset[0]==0)
+         m_faceCount[0]=m_NE[1];
+     else
+         m_faceCount[0]=0;
+     //right
+     if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
+         m_faceCount[1]=m_NE[1];
+     else
+         m_faceCount[1]=0;
+     //bottom
+     if (m_offset[1]==0)
+         m_faceCount[2]=m_NE[0];
+     else
+         m_faceCount[2]=0;
+     //top
+     if (m_mpiInfo->rank/m_NX[0]==m_NX[1]-1)
+         m_faceCount[3]=m_NE[0];
+     else
+         m_faceCount[3]=0;
      m_faceId.resize(getNumFaceElements());
+     const index_t left = (m_offset[0]==0 ? 0 : 1);
+     const index_t bottom = (m_offset[1]==0 ? 0 : 1);
+     const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
+     const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
+ #define globalNodeId(x,y) \
+     ((m_offset[0]+x)/nDOF0)*nDOF0*nDOF1+(m_offset[0]+x)%nDOF0 \
+     + ((m_offset[1]+y)/nDOF1)*nDOF0*nDOF1*m_NX[0]+((m_offset[1]+y)%nDOF1)*nDOF0
+     // set corner id's outside the parallel region
+     m_nodeId[0] = globalNodeId(0, 0);
+     m_nodeId[m_NN[0]-1] = globalNodeId(m_NN[0]-1, 0);
+     m_nodeId[m_NN[0]*(m_NN[1]-1)] = globalNodeId(0, m_NN[1]-1);
+     m_nodeId[m_NN[0]*m_NN[1]-1] = globalNodeId(m_NN[0]-1,m_NN[1]-1);
+ #undef globalNodeId
  #pragma omp parallel
      {
-         // nodes
+         // populate degrees of freedom and own nodes (identical id)
  #pragma omp for nowait
-         for (dim_t i1=0; i1<m_N1; i1++) {
+         for (dim_t i=0; i<nDOF1; i++) {
-             for (dim_t i0=0; i0<m_N0; i0++) {
+             for (dim_t j=0; j<nDOF0; j++) {
-                 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];
+                 const index_t dofIdx=j+i*nDOF0;
+                 m_dofId[dofIdx] = m_nodeId[nodeIdx]
+                     = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
              }
          }
-         // degrees of freedom
+         // populate the rest of the nodes (shared with other ranks)
+         if (m_faceCount[0]==0) { // left column
+ #pragma omp for nowait
+             for (dim_t i=0; i<nDOF1; i++) {
+                 const index_t nodeIdx=(i+bottom)*m_NN[0];
+                 const index_t dofId=(i+1)*nDOF0-1;
+                 m_nodeId[nodeIdx]
+                     = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
+             }
+         }
+         if (m_faceCount[1]==0) { // right column
+ #pragma omp for nowait
+             for (dim_t i=0; i<nDOF1; i++) {
+                 const index_t nodeIdx=(i+bottom+1)*m_NN[0]-1;
+                 const index_t dofId=i*nDOF0;
+                 m_nodeId[nodeIdx]
+                     = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
+             }
+         }
+         if (m_faceCount[2]==0) { // bottom row
+ #pragma omp for nowait
+             for (dim_t i=0; i<nDOF0; i++) {
+                 const index_t nodeIdx=i+left;
+                 const index_t dofId=nDOF0*(nDOF1-1)+i;
+                 m_nodeId[nodeIdx]
+                     = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
+             }
+         }
+         if (m_faceCount[3]==0) { // top row
  #pragma omp for nowait
-         for (dim_t k=0; k<numDOF; k++)
+             for (dim_t i=0; i<nDOF0; i++) {
-             m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;
+                 const index_t nodeIdx=m_NN[0]*(m_NN[1]-1)+i+left;
+                 const index_t dofId=i;
+                 m_nodeId[nodeIdx]
+                     = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
+             }
+         }
-         // elements
+         // populate element id's
  #pragma omp for nowait
-         for (dim_t i1=0; i1<m_NE1; i1++) {
+         for (dim_t i1=0; i1<m_NE[1]; i1++) {
-             for (dim_t i0=0; i0<m_NE0; i0++) {
+             for (dim_t i0=0; i0<m_NE[0]; i0++) {
-                 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;
              }
          }
-Line 1325 
 void Rectangle::populateSampleIds()
+Line 1590 
 void Rectangle::populateSampleIds()
      updateTagsInUse(Elements);
      // generate face offset vector and set face tags
-     const IndexVector facesPerEdge = getNumFacesPerBoundary();
      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
-     m_faceOffset.assign(facesPerEdge.size(), -1);
+     m_faceOffset.assign(4, -1);
      m_faceTags.clear();
      index_t offset=0;
-     for (size_t i=0; i<facesPerEdge.size(); i++) {
+     for (size_t i=0; i<4; i++) {
-         if (facesPerEdge[i]>0) {
+         if (m_faceCount[i]>0) {
              m_faceOffset[i]=offset;
-             offset+=facesPerEdge[i];
+             offset+=m_faceCount[i];
-             m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);
+             m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
          }
      }
      setTagMap("left", LEFT);
-Line 1348 
 void Rectangle::populateSampleIds()
+Line 1612 
 void Rectangle::populateSampleIds()
  //private
  void Rectangle::createPattern()
  {
-     const dim_t nDOF0 = (m_gNE0+1)/m_NX;
+     const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
-     const dim_t nDOF1 = (m_gNE1+1)/m_NY;
+     const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
-     const index_t left = (m_offset0==0 ? 0 : 1);
+     const index_t left = (m_offset[0]==0 ? 0 : 1);
-     const index_t bottom = (m_offset1==0 ? 0 : 1);
+     const index_t bottom = (m_offset[1]==0 ? 0 : 1);
      // populate node->DOF mapping with own degrees of freedom.
      // The rest is assigned in the loop further down
-Line 1359 
 void Rectangle::createPattern()
+Line 1623 
 void Rectangle::createPattern()
  #pragma omp parallel for
      for (index_t i=bottom; i<bottom+nDOF1; i++) {
          for (index_t j=left; j<left+nDOF0; j++) {
-             m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
+             m_dofMap[i*m_NN[0]+j]=(i-bottom)*nDOF0+j-left;
          }
      }
-Line 1371 
 void Rectangle::createPattern()
+Line 1635 
 void Rectangle::createPattern()
      RankVector neighbour;
      IndexVector offsetInShared(1,0);
      IndexVector sendShared, recvShared;
-     int numShared=0;
+     int numShared=0, expectedShared=0;
-     const int x=m_mpiInfo->rank%m_NX;
+     const int x=m_mpiInfo->rank%m_NX[0];
-     const int y=m_mpiInfo->rank/m_NX;
+     const int y=m_mpiInfo->rank/m_NX[0];
+     if (x > 0)
+         expectedShared += nDOF1;
+     if (x < m_NX[0] - 1)
+         expectedShared += nDOF1;
+     if (y > 0)
+         expectedShared += nDOF0;
+     if (y < m_NX[1] - 1)
+         expectedShared += nDOF0;
+     if (x > 0 && y > 0) expectedShared++;
+     if (x > 0 && y < m_NX[1] - 1) expectedShared++;
+     if (x < m_NX[0] - 1 && y > 0) expectedShared++;
+     if (x < m_NX[0] - 1 && y < m_NX[1] - 1) expectedShared++;
+     vector<IndexVector> rowIndices(expectedShared);
      for (int i1=-1; i1<2; i1++) {
          for (int i0=-1; i0<2; i0++) {
              // skip this rank
-Line 1382 
 void Rectangle::createPattern()
+Line 1661 
 void Rectangle::createPattern()
              // location of neighbour rank
              const int nx=x+i0;
              const int ny=y+i1;
-             if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {
+             if (nx>=0 && ny>=0 && nx<m_NX[0] && ny<m_NX[1]) {
-                 neighbour.push_back(ny*m_NX+nx);
+                 neighbour.push_back(ny*m_NX[0]+nx);
                  if (i0==0) {
                      // sharing top or bottom edge
                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
-                     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);
                      offsetInShared.push_back(offsetInShared.back()+nDOF0);
                      for (dim_t i=0; i<nDOF0; i++, numShared++) {
                          sendShared.push_back(firstDOF+i);
                          recvShared.push_back(numDOF+numShared);
                          if (i>0)
-                             colIndices[firstDOF+i-1].push_back(numShared);
+                             doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
-                         colIndices[firstDOF+i].push_back(numShared);
+                         doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
                          if (i<nDOF0-1)
-                             colIndices[firstDOF+i+1].push_back(numShared);
+                             doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
                          m_dofMap[firstNode+i]=numDOF+numShared;
                      }
                  } else if (i1==0) {
                      // sharing left or right edge
                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
-                     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);
                      offsetInShared.push_back(offsetInShared.back()+nDOF1);
                      for (dim_t i=0; i<nDOF1; i++, numShared++) {
                          sendShared.push_back(firstDOF+i*nDOF0);
                          recvShared.push_back(numDOF+numShared);
                          if (i>0)
-                             colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);
+                             doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
-                         colIndices[firstDOF+i*nDOF0].push_back(numShared);
+                         doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
                          if (i<nDOF1-1)
-                             colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
+                             doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
-                         m_dofMap[firstNode+i*m_N0]=numDOF+numShared;
+                         m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
                      }
                  } else {
                      // sharing a node
                      const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
-                     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);
                      offsetInShared.push_back(offsetInShared.back()+1);
                      sendShared.push_back(dof);
                      recvShared.push_back(numDOF+numShared);
-                     colIndices[dof].push_back(numShared);
+                     doublyLink(colIndices, rowIndices, dof, numShared);
                      m_dofMap[node]=numDOF+numShared;
                      ++numShared;
                  }
              }
          }
      }
+ #pragma omp parallel for
+     for (int i = 0; i < numShared; i++) {
+         std::sort(rowIndices[i].begin(), rowIndices[i].end());
+     }
      // create connector
      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
-Line 1443 
 void Rectangle::createPattern()
+Line 1727 
 void Rectangle::createPattern()
      // create main and couple blocks
      Paso_Pattern *mainPattern = createMainPattern();
      Paso_Pattern *colPattern, *rowPattern;
-     createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
+     createCouplePatterns(colIndices, rowIndices, numShared, &colPattern, &rowPattern);
      // allocate paso distribution
      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
-Line 1527 
 void Rectangle::addToMatrixAndRHS(Paso_S
+Line 1811 
 void Rectangle::addToMatrixAndRHS(Paso_S
      IndexVector rowIndex;
      rowIndex.push_back(m_dofMap[firstNode]);
      rowIndex.push_back(m_dofMap[firstNode+1]);
-     rowIndex.push_back(m_dofMap[firstNode+m_N0]);
+     rowIndex.push_back(m_dofMap[firstNode+m_NN[0]]);
-     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
+     rowIndex.push_back(m_dofMap[firstNode+m_NN[0]+1]);
      if (addF) {
          double *F_p=F.getSampleDataRW(0);
          for (index_t i=0; i<rowIndex.size(); i++) {
-Line 1546 
 void Rectangle::addToMatrixAndRHS(Paso_S
+Line 1830 
 void Rectangle::addToMatrixAndRHS(Paso_S
  //protected
  void Rectangle::interpolateNodesOnElements(escript::Data& out,
-                                         escript::Data& in, bool reduced) const
+                                            const escript::Data& in,
+                                            bool reduced) const
  {
      const dim_t numComp = in.getDataPointSize();
      if (reduced) {
-Line 1559 
 void Rectangle::interpolateNodesOnElemen
+Line 1844 
 void Rectangle::interpolateNodesOnElemen
              vector<double> f_10(numComp);
              vector<double> f_11(numComp);
  #pragma omp for
-             for (index_t k1=0; k1 < m_NE1; ++k1) {
+             for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
-                     double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
+                     double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
                      } /* end of component loop i */
-Line 1584 
 void Rectangle::interpolateNodesOnElemen
+Line 1869 
 void Rectangle::interpolateNodesOnElemen
              vector<double> f_10(numComp);
              vector<double> f_11(numComp);
  #pragma omp for
-             for (index_t k1=0; k1 < m_NE1; ++k1) {
+             for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
-                     double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
+                     double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
                          o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
-Line 1604 
 void Rectangle::interpolateNodesOnElemen
+Line 1889 
 void Rectangle::interpolateNodesOnElemen
  }
  //protected
- void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,
+ void Rectangle::interpolateNodesOnFaces(escript::Data& out,
+                                         const escript::Data& in,
                                          bool reduced) const
  {
      const dim_t numComp = in.getDataPointSize();
      if (reduced) {
          out.requireWrite();
-         const double c0 = 0.5;
  #pragma omp parallel
          {
              vector<double> f_00(numComp);
-Line 1619 
 void Rectangle::interpolateNodesOnFaces(
+Line 1904 
 void Rectangle::interpolateNodesOnFaces(
              vector<double> f_11(numComp);
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                 for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
                      for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
+                         o[INDEX2(i,numComp,0)] = (f_00[i] + f_01[i])/2;
                      } /* end of component loop i */
                  } /* end of k1 loop */
              } /* end of face 0 */
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                 for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
                      for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
+                         o[INDEX2(i,numComp,0)] = (f_10[i] + f_11[i])/2;
                      } /* end of component loop i */
                  } /* end of k1 loop */
              } /* end of face 1 */
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
                      for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
+                         o[INDEX2(i,numComp,0)] = (f_00[i] + f_10[i])/2;
                      } /* end of component loop i */
                  } /* end of k0 loop */
              } /* end of face 2 */
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
                      for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
+                         o[INDEX2(i,numComp,0)] = (f_01[i] + f_11[i])/2;
                      } /* end of component loop i */
                  } /* end of k0 loop */
              } /* end of face 3 */
-Line 1674 
 void Rectangle::interpolateNodesOnFaces(
+Line 1959 
 void Rectangle::interpolateNodesOnFaces(
              vector<double> f_11(numComp);
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                 for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c0*f_01[i] + c1*f_00[i];
-Line 1686 
 void Rectangle::interpolateNodesOnFaces(
+Line 1971 
 void Rectangle::interpolateNodesOnFaces(
              } /* end of face 0 */
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                 for (index_t k1=0; k1 < m_NE[1]; ++k1) {
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c1*f_10[i] + c0*f_11[i];
-Line 1698 
 void Rectangle::interpolateNodesOnFaces(
+Line 1983 
 void Rectangle::interpolateNodesOnFaces(
              } /* end of face 1 */
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c0*f_10[i] + c1*f_00[i];
-Line 1710 
 void Rectangle::interpolateNodesOnFaces(
+Line 1995 
 void Rectangle::interpolateNodesOnFaces(
              } /* end of face 2 */
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
-                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
-                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
+                     memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c0*f_11[i] + c1*f_01[i];
-Line 1724 
 void Rectangle::interpolateNodesOnFaces(
+Line 2009 
 void Rectangle::interpolateNodesOnFaces(
      }
  }
- //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;
-                         }
-                     }
-                     // add to matrix (if add_EM_S) and RHS (if add_EM_F)
+ namespace
-                     const index_t firstNode=m_N0*k1+k0;
+ {
-                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
+     // Calculates a guassian blur colvolution matrix for 2D
-                 } // end k0 loop
+     // See wiki article on the subject
-             } // end k1 loop
+     double* get2DGauss(unsigned radius, double sigma)
-         } // end of colouring
+     {
-     } // end of parallel region
+         double* arr=new double[(radius*2+1)*(radius*2+1)];
+         double common=M_1_PI*0.5*1/(sigma*sigma);
+         double total=0;
+         int r=static_cast<int>(radius);
+         for (int y=-r;y<=r;++y)
+         {
+             for (int x=-r;x<=r;++x)
+             {
+                 arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));
+                 total+=arr[(x+r)+(y+r)*(r*2+1)];
+             }
+         }
+         double invtotal=1/total;
+         for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)
+         {
+             arr[p]*=invtotal;
+         }
+         return arr;
+     }
+     // applies conv to source to get a point.
+     // (xp, yp) are the coords in the source matrix not the destination matrix
+     double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)
+     {
+         size_t bx=xp-radius, by=yp-radius;
+         size_t sbase=bx+by*width;
+         double result=0;
+         for (int y=0;y<2*radius+1;++y)
+         {
+             for (int x=0;x<2*radius+1;++x)
+             {
+                 result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];
+             }
+         }
+         return result;
+     }
  }
- //protected
- 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
- {
-     const double h0 = m_l0/m_gNE0;
-     const double h1 = m_l1/m_gNE1;
-     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();
+ /* This is a wrapper for filtered (and non-filtered) randoms
- #pragma omp parallel
+  * For detailed doco see randomFillWorker
+ */
+ escript::Data Rectangle::randomFill(const escript::DataTypes::ShapeType& shape,
+        const escript::FunctionSpace& what,
+        long seed, const boost::python::tuple& filter) const
+ {
+     int numvals=escript::DataTypes::noValues(shape);
+     if (len(filter)>0 && (numvals!=1))
      {
-         for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
+         throw RipleyException("Ripley only supports filters for scalar data.");
- #pragma omp for
+     }
-             for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
+     escript::Data res=randomFillWorker(shape, seed, filter);
-                 for (index_t k0=0; k0<m_NE0; ++k0)  {
+     if (res.getFunctionSpace()!=what)
-                     bool add_EM_S=false;
+     {
-                     bool add_EM_F=false;
+         escript::Data r=escript::Data(res, what);
-                     vector<double> EM_S(4*4, 0);
+         return r;
-                     vector<double> EM_F(4, 0);
+     }
-                     const index_t e = k0 + m_NE0*k1;
+     return res;
-                     ///////////////
-                     // process A //
-                     ///////////////
-                     if (!A.isEmpty()) {
-                         add_EM_S=true;
-                         const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
-                         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_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[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
-                         EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
-                         EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
-                         EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
-                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
-                         EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
-                         EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
-                         EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
-                         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
  }
- //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;
-     rhs.requireWrite();
+ /* This routine produces a Data object filled with smoothed random data.
- #pragma omp parallel
+ The dimensions of the rectangle being filled are internal[0] x internal[1] points.
-     {
+ A parameter radius  gives the size of the stencil used for the smoothing.
-         for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
+ A point on the left hand edge for example, will still require `radius` extra points to the left
- #pragma omp for
+ in order to complete the stencil.
-             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;
-                             }
-                         }
-                     }
-                     // add to matrix (if add_EM_S) and RHS (if add_EM_F)
+ All local calculation is done on an array called `src`, with
-                     const index_t firstNode=m_N0*k1+k0;
+ dimensions = ext[0] * ext[1].
-                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
+ Where ext[i]= internal[i]+2*radius.
-                             firstNode, numEq, numComp);
-                 } // end k0 loop
-             } // end k1 loop
-         } // end of colouring
-     } // end of parallel region
- }
- //protected
+ Now for MPI there is overlap to deal with. We need to share both the overlapping
- void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,
+ values themselves but also the external region.
-         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 = -.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();
+ In a hypothetical 1-D case:
- #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);
-                         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;