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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4370 by caltinay, Fri Apr 19 06:15:24 2013 UTC revision 4705 by jfenwick, Fri Feb 21 02:36:15 2014 UTC
# Line 1  Line 1 
1    
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2013 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
# Line 9  Line 9 
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  * Development since 2012 by School of Earth Sciences  * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
17  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
18  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
19  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
20    #include <ripley/DefaultAssembler2D.h>
21    #include <ripley/WaveAssembler2D.h>
22    #include <boost/scoped_array.hpp>
23    #include "esysUtils/EsysRandom.h"
24    #include "blocktools.h"
25    
26  #ifdef USE_NETCDF  #ifdef USE_NETCDF
27  #include <netcdfcpp.h>  #include <netcdfcpp.h>
# Line 36  using esysUtils::FileWriter; Line 42  using esysUtils::FileWriter;
42  namespace ripley {  namespace ripley {
43    
44  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
45                       double y1, int d0, int d1) :                       double y1, int d0, int d1,
46                         const std::vector<double>& points,
47                         const std::vector<int>& tags,
48                         const simap_t& tagnamestonums) :
49      RipleyDomain(2)      RipleyDomain(2)
50  {  {
51      // ignore subdivision parameters for serial run      // ignore subdivision parameters for serial run
# Line 138  Rectangle::Rectangle(int n0, int n1, dou Line 147  Rectangle::Rectangle(int n0, int n1, dou
147    
148      populateSampleIds();      populateSampleIds();
149      createPattern();      createPattern();
150        assembler = new DefaultAssembler2D(this, m_dx, m_NX, m_NE, m_NN);
151        for (map<string, int>::const_iterator i = tagnamestonums.begin();
152                i != tagnamestonums.end(); i++) {
153            setTagMap(i->first, i->second);
154        }
155        addPoints(tags.size(), &points[0], &tags[0]);
156  }  }
157    
158  Rectangle::~Rectangle()  Rectangle::~Rectangle()
159  {  {
160      Paso_SystemMatrixPattern_free(m_pattern);      Paso_SystemMatrixPattern_free(m_pattern);
161      Paso_Connector_free(m_connector);      Paso_Connector_free(m_connector);
162        delete assembler;
163  }  }
164    
165  string Rectangle::getDescription() const  string Rectangle::getDescription() const
# Line 166  bool Rectangle::operator==(const Abstrac Line 182  bool Rectangle::operator==(const Abstrac
182  }  }
183    
184  void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,  void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
185              const vector<int>& first, const vector<int>& numValues,              const ReaderParameters& params) const
             const vector<int>& multiplier) const  
186  {  {
187  #ifdef USE_NETCDF  #ifdef USE_NETCDF
188      // check destination function space      // check destination function space
# Line 182  void Rectangle::readNcGrid(escript::Data Line 197  void Rectangle::readNcGrid(escript::Data
197      } else      } else
198          throw RipleyException("readNcGrid(): invalid function space for output data object");          throw RipleyException("readNcGrid(): invalid function space for output data object");
199    
200      if (first.size() != 2)      if (params.first.size() != 2)
201          throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");          throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
202    
203      if (numValues.size() != 2)      if (params.numValues.size() != 2)
204          throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");          throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
205    
206      if (multiplier.size() != 2)      if (params.multiplier.size() != 2)
207          throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");          throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");
208      for (size_t i=0; i<multiplier.size(); i++)      for (size_t i=0; i<params.multiplier.size(); i++)
209          if (multiplier[i]<1)          if (params.multiplier[i]<1)
210              throw RipleyException("readNcGrid(): all multipliers must be positive");              throw RipleyException("readNcGrid(): all multipliers must be positive");
211        if (params.reverse.size() != 2)
212            throw RipleyException("readNcGrid(): argument 'reverse' must have 2 entries");
213    
214      // check file existence and size      // check file existence and size
215      NcFile f(filename.c_str(), NcFile::ReadOnly);      NcFile f(filename.c_str(), NcFile::ReadOnly);
# Line 209  void Rectangle::readNcGrid(escript::Data Line 226  void Rectangle::readNcGrid(escript::Data
226          throw RipleyException("readNcGrid(): only scalar data supported");          throw RipleyException("readNcGrid(): only scalar data supported");
227    
228      const int dims = var->num_dims();      const int dims = var->num_dims();
229      const long *edges = var->edges();      boost::scoped_array<long> edges(var->edges());
230    
231      // is this a slice of the data object (dims!=2)?      // is this a slice of the data object (dims!=2)?
232      // note the expected ordering of edges (as in numpy: y,x)      // note the expected ordering of edges (as in numpy: y,x)
233      if ( (dims==2 && (numValues[1] > edges[0] || numValues[0] > edges[1]))      if ( (dims==2 && (params.numValues[1] > edges[0] || params.numValues[0] > edges[1]))
234              || (dims==1 && numValues[1]>1) ) {              || (dims==1 && params.numValues[1]>1) ) {
235          throw RipleyException("readNcGrid(): not enough data in file");          throw RipleyException("readNcGrid(): not enough data in file");
236      }      }
237    
238      // check if this rank contributes anything      // check if this rank contributes anything
239      if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset[0] ||      if (params.first[0] >= m_offset[0]+myN0 ||
240              first[1] >= m_offset[1]+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset[1])              params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
241                params.first[1] >= m_offset[1]+myN1 ||
242                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1])
243          return;          return;
244    
245      // now determine how much this rank has to write      // now determine how much this rank has to write
246    
247      // first coordinates in data object to write to      // first coordinates in data object to write to
248      const int first0 = max(0, first[0]-m_offset[0]);      const int first0 = max(0, params.first[0]-m_offset[0]);
249      const int first1 = max(0, first[1]-m_offset[1]);      const int first1 = max(0, params.first[1]-m_offset[1]);
250      // indices to first value in file      // indices to first value in file (not accounting for reverse yet)
251      const int idx0 = max(0, m_offset[0]-first[0]);      int idx0 = max(0, m_offset[0]-params.first[0]);
252      const int idx1 = max(0, m_offset[1]-first[1]);      int idx1 = max(0, m_offset[1]-params.first[1]);
253      // number of values to read      // number of values to read
254      const int num0 = min(numValues[0]-idx0, myN0-first0);      const int num0 = min(params.numValues[0]-idx0, myN0-first0);
255      const int num1 = min(numValues[1]-idx1, myN1-first1);      const int num1 = min(params.numValues[1]-idx1, myN1-first1);
256    
257        // make sure we read the right block if going backwards through file
258        if (params.reverse[0])
259            idx0 = edges[dims-1]-num0-idx0;
260        if (dims>1 && params.reverse[1])
261            idx1 = edges[dims-2]-num1-idx1;
262    
263      vector<double> values(num0*num1);      vector<double> values(num0*num1);
264      if (dims==2) {      if (dims==2) {
# Line 247  void Rectangle::readNcGrid(escript::Data Line 272  void Rectangle::readNcGrid(escript::Data
272      const int dpp = out.getNumDataPointsPerSample();      const int dpp = out.getNumDataPointsPerSample();
273      out.requireWrite();      out.requireWrite();
274    
275        // helpers for reversing
276        const int x0 = (params.reverse[0] ? num0-1 : 0);
277        const int x_mult = (params.reverse[0] ? -1 : 1);
278        const int y0 = (params.reverse[1] ? num1-1 : 0);
279        const int y_mult = (params.reverse[1] ? -1 : 1);
280    
281      for (index_t y=0; y<num1; y++) {      for (index_t y=0; y<num1; y++) {
282  #pragma omp parallel for  #pragma omp parallel for
283          for (index_t x=0; x<num0; x++) {          for (index_t x=0; x<num0; x++) {
284              const int baseIndex = first0+x*multiplier[0]              const int baseIndex = first0+x*params.multiplier[0]
285                                    +(first1+y*multiplier[1])*myN0;                                    +(first1+y*params.multiplier[1])*myN0;
286              const int srcIndex = y*num0+x;              const int srcIndex = (y0+y_mult*y)*num0+(x0+x_mult*x);
287              if (!isnan(values[srcIndex])) {              if (!isnan(values[srcIndex])) {
288                  for (index_t m1=0; m1<multiplier[1]; m1++) {                  for (index_t m1=0; m1<params.multiplier[1]; m1++) {
289                      for (index_t m0=0; m0<multiplier[0]; m0++) {                      for (index_t m0=0; m0<params.multiplier[0]; m0++) {
290                          const int dataIndex = baseIndex+m0+m1*myN0;                          const int dataIndex = baseIndex+m0+m1*myN0;
291                          double* dest = out.getSampleDataRW(dataIndex);                          double* dest = out.getSampleDataRW(dataIndex);
292                          for (index_t q=0; q<dpp; q++) {                          for (index_t q=0; q<dpp; q++) {
# Line 272  void Rectangle::readNcGrid(escript::Data Line 303  void Rectangle::readNcGrid(escript::Data
303  }  }
304    
305  void Rectangle::readBinaryGrid(escript::Data& out, string filename,  void Rectangle::readBinaryGrid(escript::Data& out, string filename,
306                                 const vector<int>& first,                                 const ReaderParameters& params) const
307                                 const vector<int>& numValues,  {
308                                 const vector<int>& multiplier) const      // the mapping is not universally correct but should work on our
309        // supported platforms
310        switch (params.dataType) {
311            case DATATYPE_INT32:
312                readBinaryGridImpl<int>(out, filename, params);
313                break;
314            case DATATYPE_FLOAT32:
315                readBinaryGridImpl<float>(out, filename, params);
316                break;
317            case DATATYPE_FLOAT64:
318                readBinaryGridImpl<double>(out, filename, params);
319                break;
320            default:
321                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
322        }
323    }
324    
325    template<typename ValueType>
326    void Rectangle::readBinaryGridImpl(escript::Data& out, const string& filename,
327                                       const ReaderParameters& params) const
328  {  {
329      // check destination function space      // check destination function space
330      int myN0, myN1;      int myN0, myN1;
# Line 296  void Rectangle::readBinaryGrid(escript:: Line 346  void Rectangle::readBinaryGrid(escript::
346      f.seekg(0, ios::end);      f.seekg(0, ios::end);
347      const int numComp = out.getDataPointSize();      const int numComp = out.getDataPointSize();
348      const int filesize = f.tellg();      const int filesize = f.tellg();
349      const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(float);      const int reqsize = params.numValues[0]*params.numValues[1]*numComp*sizeof(ValueType);
350      if (filesize < reqsize) {      if (filesize < reqsize) {
351          f.close();          f.close();
352          throw RipleyException("readBinaryGrid(): not enough data in file");          throw RipleyException("readBinaryGrid(): not enough data in file");
353      }      }
354    
355      // check if this rank contributes anything      // check if this rank contributes anything
356      if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0] <= m_offset[0] ||      if (params.first[0] >= m_offset[0]+myN0 ||
357              first[1] >= m_offset[1]+myN1 || first[1]+numValues[1] <= m_offset[1]) {              params.first[0]+params.numValues[0] <= m_offset[0] ||
358                params.first[1] >= m_offset[1]+myN1 ||
359                params.first[1]+params.numValues[1] <= m_offset[1]) {
360          f.close();          f.close();
361          return;          return;
362      }      }
# Line 312  void Rectangle::readBinaryGrid(escript:: Line 364  void Rectangle::readBinaryGrid(escript::
364      // now determine how much this rank has to write      // now determine how much this rank has to write
365    
366      // first coordinates in data object to write to      // first coordinates in data object to write to
367      const int first0 = max(0, first[0]-m_offset[0]);      const int first0 = max(0, params.first[0]-m_offset[0]);
368      const int first1 = max(0, first[1]-m_offset[1]);      const int first1 = max(0, params.first[1]-m_offset[1]);
369      // indices to first value in file      // indices to first value in file
370      const int idx0 = max(0, m_offset[0]-first[0]);      const int idx0 = max(0, m_offset[0]-params.first[0]);
371      const int idx1 = max(0, m_offset[1]-first[1]);      const int idx1 = max(0, m_offset[1]-params.first[1]);
372      // number of values to read      // number of values to read
373      const int num0 = min(numValues[0]-idx0, myN0-first0);      const int num0 = min(params.numValues[0]-idx0, myN0-first0);
374      const int num1 = min(numValues[1]-idx1, myN1-first1);      const int num1 = min(params.numValues[1]-idx1, myN1-first1);
375    
376      out.requireWrite();      out.requireWrite();
377      vector<float> values(num0*numComp);      vector<ValueType> values(num0*numComp);
378      const int dpp = out.getNumDataPointsPerSample();      const int dpp = out.getNumDataPointsPerSample();
379    
380      for (index_t y=0; y<num1; y++) {      for (int y=0; y<num1; y++) {
381          const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);          const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]);
382          f.seekg(fileofs*sizeof(float));          f.seekg(fileofs*sizeof(ValueType));
383          f.read((char*)&values[0], num0*numComp*sizeof(float));          f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
384          for (index_t x=0; x<num0; x++) {          for (int x=0; x<num0; x++) {
385              const int baseIndex = first0+x*multiplier[0]              const int baseIndex = first0+x*params.multiplier[0]
386                                      +(first1+y*multiplier[1])*myN0;                                      +(first1+y*params.multiplier[1])*myN0;
387              for (index_t m1=0; m1<multiplier[1]; m1++) {              for (int m1=0; m1<params.multiplier[1]; m1++) {
388                  for (index_t m0=0; m0<multiplier[0]; m0++) {                  for (int m0=0; m0<params.multiplier[0]; m0++) {
389                      const int dataIndex = baseIndex+m0+m1*myN0;                      const int dataIndex = baseIndex+m0+m1*myN0;
390                      double* dest = out.getSampleDataRW(dataIndex);                      double* dest = out.getSampleDataRW(dataIndex);
391                      for (index_t c=0; c<numComp; c++) {                      for (int c=0; c<numComp; c++) {
392                          if (!std::isnan(values[x*numComp+c])) {                          ValueType val = values[x*numComp+c];
393                              for (index_t q=0; q<dpp; q++) {  
394                                  *dest++ = static_cast<double>(values[x*numComp+c]);                          if (params.byteOrder != BYTEORDER_NATIVE) {
395                                char* cval = reinterpret_cast<char*>(&val);
396                                // this will alter val!!
397                                byte_swap32(cval);
398                            }
399                            if (!std::isnan(val)) {
400                                for (int q=0; q<dpp; q++) {
401                                    *dest++ = static_cast<double>(val);
402                              }                              }
403                          }                          }
404                      }                      }
# Line 398  void Rectangle::writeBinaryGridImpl(cons Line 457  void Rectangle::writeBinaryGridImpl(cons
457      if (numComp > 1 || dpp > 1)      if (numComp > 1 || dpp > 1)
458          throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");          throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
459    
     escript::Data* _in = const_cast<escript::Data*>(&in);  
460      const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;      const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;
461    
462      // from here on we know that each sample consists of one value      // from here on we know that each sample consists of one value
463      FileWriter* fw = new FileWriter();      FileWriter fw;
464      fw->openFile(filename, fileSize);      fw.openFile(filename, fileSize);
465      MPIBarrier();      MPIBarrier();
466    
467      for (index_t y=0; y<myN1; y++) {      for (index_t y=0; y<myN1; y++) {
# Line 411  void Rectangle::writeBinaryGridImpl(cons Line 469  void Rectangle::writeBinaryGridImpl(cons
469          ostringstream oss;          ostringstream oss;
470    
471          for (index_t x=0; x<myN0; x++) {          for (index_t x=0; x<myN0; x++) {
472              const double* sample = _in->getSampleDataRO(y*myN0+x);              const double* sample = in.getSampleDataRO(y*myN0+x);
473              ValueType fvalue = static_cast<ValueType>(*sample);              ValueType fvalue = static_cast<ValueType>(*sample);
474              if (byteOrder == BYTEORDER_NATIVE) {              if (byteOrder == BYTEORDER_NATIVE) {
475                  oss.write((char*)&fvalue, sizeof(fvalue));                  oss.write((char*)&fvalue, sizeof(fvalue));
# Line 420  void Rectangle::writeBinaryGridImpl(cons Line 478  void Rectangle::writeBinaryGridImpl(cons
478                  oss.write(byte_swap32(value), sizeof(fvalue));                  oss.write(byte_swap32(value), sizeof(fvalue));
479              }              }
480          }          }
481          fw->writeAt(oss, fileofs);          fw.writeAt(oss, fileofs);
482      }      }
483      fw->close();      fw.close();
484  }  }
485    
486  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
# Line 582  const int* Rectangle::borrowSampleRefere Line 640  const int* Rectangle::borrowSampleRefere
640          case FaceElements:          case FaceElements:
641          case ReducedFaceElements:          case ReducedFaceElements:
642              return &m_faceId[0];              return &m_faceId[0];
643            case Points:
644                return &m_diracPointNodeIDs[0];
645          default:          default:
646              break;              break;
647      }      }
# Line 839  void Rectangle::assembleCoordinates(escr Line 899  void Rectangle::assembleCoordinates(escr
899  }  }
900    
901  //protected  //protected
902  void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const  void Rectangle::assembleGradient(escript::Data& out, const escript::Data& in) const
903  {  {
904      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
905      const double cx0 = -1./m_dx[0];      const double cx0 = .21132486540518711775/m_dx[0];
906      const double cx1 = -.78867513459481288225/m_dx[0];      const double cx1 = .78867513459481288225/m_dx[0];
907      const double cx2 = -.5/m_dx[0];      const double cx2 = 1./m_dx[0];
908      const double cx3 = -.21132486540518711775/m_dx[0];      const double cy0 = .21132486540518711775/m_dx[1];
909      const double cx4 = .21132486540518711775/m_dx[0];      const double cy1 = .78867513459481288225/m_dx[1];
910      const double cx5 = .5/m_dx[0];      const double cy2 = 1./m_dx[1];
     const double cx6 = .78867513459481288225/m_dx[0];  
     const double cx7 = 1./m_dx[0];  
     const double cy0 = -1./m_dx[1];  
     const double cy1 = -.78867513459481288225/m_dx[1];  
     const double cy2 = -.5/m_dx[1];  
     const double cy3 = -.21132486540518711775/m_dx[1];  
     const double cy4 = .21132486540518711775/m_dx[1];  
     const double cy5 = .5/m_dx[1];  
     const double cy6 = .78867513459481288225/m_dx[1];  
     const double cy7 = 1./m_dx[1];  
911    
912      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
913          out.requireWrite();          out.requireWrite();
# Line 876  void Rectangle::assembleGradient(escript Line 926  void Rectangle::assembleGradient(escript
926                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
927                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
928                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
929                          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;
930                          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;
931                          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;
932                          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;
933                          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;
934                          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;
935                          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;
936                          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;
937                      } // end of component loop i                      } // end of component loop i
938                  } // end of k0 loop                  } // end of k0 loop
939              } // end of k1 loop              } // end of k1 loop
# Line 905  void Rectangle::assembleGradient(escript Line 955  void Rectangle::assembleGradient(escript
955                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
956                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
957                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
958                          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;
959                          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;
960                      } // end of component loop i                      } // end of component loop i
961                  } // end of k0 loop                  } // end of k0 loop
962              } // end of k1 loop              } // end of k1 loop
# Line 928  void Rectangle::assembleGradient(escript Line 978  void Rectangle::assembleGradient(escript
978                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
979                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
980                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
981                          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;
982                          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;
983                          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;
984                          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;
985                      } // end of component loop i                      } // end of component loop i
986                  } // end of k1 loop                  } // end of k1 loop
987              } // end of face 0              } // end of face 0
# Line 944  void Rectangle::assembleGradient(escript Line 994  void Rectangle::assembleGradient(escript
994                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
995                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
996                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
997                          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;
998                          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;
999                          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;
1000                          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;
1001                      } // end of component loop i                      } // end of component loop i
1002                  } // end of k1 loop                  } // end of k1 loop
1003              } // end of face 1              } // end of face 1
# Line 960  void Rectangle::assembleGradient(escript Line 1010  void Rectangle::assembleGradient(escript
1010                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
1011                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1012                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1013                          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;
1014                          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;
1015                          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;
1016                          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;
1017                      } // end of component loop i                      } // end of component loop i
1018                  } // end of k0 loop                  } // end of k0 loop
1019              } // end of face 2              } // end of face 2
# Line 976  void Rectangle::assembleGradient(escript Line 1026  void Rectangle::assembleGradient(escript
1026                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1027                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1028                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1029                          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;
1030                          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;
1031                          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;
1032                          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;
1033                      } // end of component loop i                      } // end of component loop i
1034                  } // end of k0 loop                  } // end of k0 loop
1035              } // end of face 3              } // end of face 3
# Line 1002  void Rectangle::assembleGradient(escript Line 1052  void Rectangle::assembleGradient(escript
1052                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
1053                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1054                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1055                          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;
1056                          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;
1057                      } // end of component loop i                      } // end of component loop i
1058                  } // end of k1 loop                  } // end of k1 loop
1059              } // end of face 0              } // end of face 0
# Line 1016  void Rectangle::assembleGradient(escript Line 1066  void Rectangle::assembleGradient(escript
1066                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1067                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1068                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1069                          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;
1070                          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;
1071                      } // end of component loop i                      } // end of component loop i
1072                  } // end of k1 loop                  } // end of k1 loop
1073              } // end of face 1              } // end of face 1
# Line 1030  void Rectangle::assembleGradient(escript Line 1080  void Rectangle::assembleGradient(escript
1080                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
1081                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1082                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1083                          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;
1084                          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;
1085                      } // end of component loop i                      } // end of component loop i
1086                  } // end of k0 loop                  } // end of k0 loop
1087              } // end of face 2              } // end of face 2
# Line 1044  void Rectangle::assembleGradient(escript Line 1094  void Rectangle::assembleGradient(escript
1094                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1095                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1096                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1097                          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;
1098                          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;
1099                      } // end of component loop i                      } // end of component loop i
1100                  } // end of k0 loop                  } // end of k0 loop
1101              } // end of face 3              } // end of face 3
# Line 1054  void Rectangle::assembleGradient(escript Line 1104  void Rectangle::assembleGradient(escript
1104  }  }
1105    
1106  //protected  //protected
1107  void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const  void Rectangle::assembleIntegrate(vector<double>& integrals,
1108                                      const escript::Data& arg) const
1109  {  {
1110      const dim_t numComp = arg.getDataPointSize();      const dim_t numComp = arg.getDataPointSize();
1111      const index_t left = (m_offset[0]==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
# Line 1240  dim_t Rectangle::insertNeighbourNodes(In Line 1291  dim_t Rectangle::insertNeighbourNodes(In
1291  }  }
1292    
1293  //protected  //protected
1294  void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const  void Rectangle::nodesToDOF(escript::Data& out, const escript::Data& in) const
1295  {  {
1296      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1297      out.requireWrite();      out.requireWrite();
# Line 1260  void Rectangle::nodesToDOF(escript::Data Line 1311  void Rectangle::nodesToDOF(escript::Data
1311  }  }
1312    
1313  //protected  //protected
1314  void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const  void Rectangle::dofToNodes(escript::Data& out, const escript::Data& in) const
1315  {  {
1316      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1317      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
1318      in.requireWrite();      // expand data object if necessary to be able to grab the whole data
1319      Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));      const_cast<escript::Data*>(&in)->expand();
1320        Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));
1321    
1322      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
1323      out.requireWrite();      out.requireWrite();
# Line 1635  void Rectangle::addToMatrixAndRHS(Paso_S Line 1687  void Rectangle::addToMatrixAndRHS(Paso_S
1687    
1688  //protected  //protected
1689  void Rectangle::interpolateNodesOnElements(escript::Data& out,  void Rectangle::interpolateNodesOnElements(escript::Data& out,
1690                                          escript::Data& in, bool reduced) const                                             const escript::Data& in,
1691                                               bool reduced) const
1692  {  {
1693      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1694      if (reduced) {      if (reduced) {
# Line 1693  void Rectangle::interpolateNodesOnElemen Line 1746  void Rectangle::interpolateNodesOnElemen
1746  }  }
1747    
1748  //protected  //protected
1749  void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,  void Rectangle::interpolateNodesOnFaces(escript::Data& out,
1750                                            const escript::Data& in,
1751                                          bool reduced) const                                          bool reduced) const
1752  {  {
1753      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1754      if (reduced) {      if (reduced) {
1755          out.requireWrite();          out.requireWrite();
         const double c0 = 0.5;  
1756  #pragma omp parallel  #pragma omp parallel
1757          {          {
1758              vector<double> f_00(numComp);              vector<double> f_00(numComp);
# Line 1713  void Rectangle::interpolateNodesOnFaces( Line 1766  void Rectangle::interpolateNodesOnFaces(
1766                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1767                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1768                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1769                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = (f_00[i] + f_01[i])/2;
1770                      } /* end of component loop i */                      } /* end of component loop i */
1771                  } /* end of k1 loop */                  } /* end of k1 loop */
1772              } /* end of face 0 */              } /* end of face 0 */
# Line 1724  void Rectangle::interpolateNodesOnFaces( Line 1777  void Rectangle::interpolateNodesOnFaces(
1777                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1778                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1779                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1780                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = (f_10[i] + f_11[i])/2;
1781                      } /* end of component loop i */                      } /* end of component loop i */
1782                  } /* end of k1 loop */                  } /* end of k1 loop */
1783              } /* end of face 1 */              } /* end of face 1 */
# Line 1735  void Rectangle::interpolateNodesOnFaces( Line 1788  void Rectangle::interpolateNodesOnFaces(
1788                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1789                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1790                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1791                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = (f_00[i] + f_10[i])/2;
1792                      } /* end of component loop i */                      } /* end of component loop i */
1793                  } /* end of k0 loop */                  } /* end of k0 loop */
1794              } /* end of face 2 */              } /* end of face 2 */
# Line 1746  void Rectangle::interpolateNodesOnFaces( Line 1799  void Rectangle::interpolateNodesOnFaces(
1799                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1800                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1801                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1802                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = (f_01[i] + f_11[i])/2;
1803                      } /* end of component loop i */                      } /* end of component loop i */
1804                  } /* end of k0 loop */                  } /* end of k0 loop */
1805              } /* end of face 3 */              } /* end of face 3 */
# Line 1813  void Rectangle::interpolateNodesOnFaces( Line 1866  void Rectangle::interpolateNodesOnFaces(
1866      }      }
1867  }  }
1868    
 //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  
 {  
     /* GENERATOR SNIP_PDE_SINGLE_PRE TOP */  
     const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];  
     const double w1 = 0.041666666666666666667;  
     const double w2 = -0.15550211698203655390;  
     const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];  
     const double w4 = -0.01116454968463011277*m_dx[1]/m_dx[0];  
     const double w5 = 0.011164549684630112770;  
     const double w6 = -0.041666666666666666667*m_dx[1]/m_dx[0];  
     const double w7 = 0.1555021169820365539*m_dx[0]/m_dx[1];  
     const double w8 = 0.01116454968463011277*m_dx[0]/m_dx[1];  
     const double w9 = -0.25000000000000000000;  
     const double w10 = -0.16666666666666666667*m_dx[0]/m_dx[1];  
     const double w11 = -0.032861463941450536761*m_dx[1];  
     const double w12 = -0.032861463941450536761*m_dx[0];  
     const double w13 = -0.12264065304058601714*m_dx[1];  
     const double w14 = -0.0023593469594139828636*m_dx[1];  
     const double w15 = -0.008805202725216129906*m_dx[0];  
     const double w16 = -0.008805202725216129906*m_dx[1];  
     const double w17 = -0.12264065304058601714*m_dx[0];  
     const double w18 = -0.0023593469594139828636*m_dx[0];  
     const double w19 = -0.16666666666666666667*m_dx[1];  
     const double w20 = -0.083333333333333333333*m_dx[0];  
     const double w21 = 0.025917019497006092316*m_dx[0]*m_dx[1];  
     const double w22 = 0.0018607582807716854616*m_dx[0]*m_dx[1];  
     const double w23 = 0.0069444444444444444444*m_dx[0]*m_dx[1];  
     const double w24 = 0.09672363354357992482*m_dx[0]*m_dx[1];  
     const double w25 = 0.000498588678642297402*m_dx[0]*m_dx[1];  
     const double w26 = 0.19716878364870322056*m_dx[1];  
     const double w27 = 0.052831216351296779436*m_dx[1];  
     const double w28 = 0.19716878364870322056*m_dx[0];  
     const double w29 = 0.052831216351296779436*m_dx[0];  
     /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */  
1869    
     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_NE[1]; k1+=2) {  
                 for (index_t k0=0; k0<m_NE[0]; ++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_NE[0]*k1;  
                     /* GENERATOR SNIP_PDE_SINGLE TOP */  
                     ///////////////  
                     // 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_01_0 = A_p[INDEX3(0,1,0,2,2)];  
                             const double A_10_0 = A_p[INDEX3(1,0,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_01_1 = A_p[INDEX3(0,1,1,2,2)];  
                             const double A_10_1 = A_p[INDEX3(1,0,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_01_2 = A_p[INDEX3(0,1,2,2,2)];  
                             const double A_10_2 = A_p[INDEX3(1,0,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_01_3 = A_p[INDEX3(0,1,3,2,2)];  
                             const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];  
                             const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];  
                             const double tmp0 = w3*(A_11_0 + A_11_1 + A_11_2 + A_11_3);  
                             const double tmp1 = w1*(A_01_0 + A_01_3 - A_10_1 - A_10_2);  
                             const double tmp2 = w4*(A_00_2 + A_00_3);  
                             const double tmp3 = w0*(A_00_0 + A_00_1);  
                             const double tmp4 = w5*(A_01_2 - A_10_3);  
                             const double tmp5 = w2*(-A_01_1 + A_10_0);  
                             const double tmp6 = w5*(A_01_3 + A_10_0);  
                             const double tmp7 = w3*(-A_11_0 - A_11_1 - A_11_2 - A_11_3);  
                             const double tmp8 = w6*(A_00_0 + A_00_1 + A_00_2 + A_00_3);  
                             const double tmp9 = w1*(A_01_1 + A_01_2 + A_10_1 + A_10_2);  
                             const double tmp10 = w2*(-A_01_0 - A_10_3);  
                             const double tmp11 = w4*(A_00_0 + A_00_1);  
                             const double tmp12 = w0*(A_00_2 + A_00_3);  
                             const double tmp13 = w5*(A_01_1 - A_10_0);  
                             const double tmp14 = w2*(-A_01_2 + A_10_3);  
                             const double tmp15 = w7*(A_11_0 + A_11_2);  
                             const double tmp16 = w4*(-A_00_2 - A_00_3);  
                             const double tmp17 = w0*(-A_00_0 - A_00_1);  
                             const double tmp18 = w5*(A_01_3 + A_10_3);  
                             const double tmp19 = w8*(A_11_1 + A_11_3);  
                             const double tmp20 = w2*(-A_01_0 - A_10_0);  
                             const double tmp21 = w7*(A_11_1 + A_11_3);  
                             const double tmp22 = w4*(-A_00_0 - A_00_1);  
                             const double tmp23 = w0*(-A_00_2 - A_00_3);  
                             const double tmp24 = w5*(A_01_0 + A_10_0);  
                             const double tmp25 = w8*(A_11_0 + A_11_2);  
                             const double tmp26 = w2*(-A_01_3 - A_10_3);  
                             const double tmp27 = w5*(-A_01_1 - A_10_2);  
                             const double tmp28 = w1*(-A_01_0 - A_01_3 - A_10_0 - A_10_3);  
                             const double tmp29 = w2*(A_01_2 + A_10_1);  
                             const double tmp30 = w7*(-A_11_1 - A_11_3);  
                             const double tmp31 = w1*(-A_01_1 - A_01_2 + A_10_0 + A_10_3);  
                             const double tmp32 = w5*(-A_01_0 + A_10_2);  
                             const double tmp33 = w8*(-A_11_0 - A_11_2);  
                             const double tmp34 = w6*(-A_00_0 - A_00_1 - A_00_2 - A_00_3);  
                             const double tmp35 = w2*(A_01_3 - A_10_1);  
                             const double tmp36 = w5*(A_01_0 + A_10_3);  
                             const double tmp37 = w2*(-A_01_3 - A_10_0);  
                             const double tmp38 = w7*(-A_11_0 - A_11_2);  
                             const double tmp39 = w5*(-A_01_3 + A_10_1);  
                             const double tmp40 = w8*(-A_11_1 - A_11_3);  
                             const double tmp41 = w2*(A_01_0 - A_10_2);  
                             const double tmp42 = w5*(A_01_1 - A_10_3);  
                             const double tmp43 = w2*(-A_01_2 + A_10_0);  
                             const double tmp44 = w5*(A_01_2 - A_10_0);  
                             const double tmp45 = w2*(-A_01_1 + A_10_3);  
                             const double tmp46 = w5*(-A_01_0 + A_10_1);  
                             const double tmp47 = w2*(A_01_3 - A_10_2);  
                             const double tmp48 = w5*(-A_01_1 - A_10_1);  
                             const double tmp49 = w2*(A_01_2 + A_10_2);  
                             const double tmp50 = w5*(-A_01_3 + A_10_2);  
                             const double tmp51 = w2*(A_01_0 - A_10_1);  
                             const double tmp52 = w5*(-A_01_2 - A_10_1);  
                             const double tmp53 = w2*(A_01_1 + A_10_2);  
                             const double tmp54 = w5*(-A_01_2 - A_10_2);  
                             const double tmp55 = w2*(A_01_1 + A_10_1);  
                             EM_S[INDEX2(0,0,4)]+=tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp9;  
                             EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;  
                             EM_S[INDEX2(0,2,4)]+=tmp31 + tmp34 + tmp38 + tmp39 + tmp40 + tmp41;  
                             EM_S[INDEX2(0,3,4)]+=tmp28 + tmp52 + tmp53 + tmp7 + tmp8;  
                             EM_S[INDEX2(1,0,4)]+=tmp0 + tmp2 + tmp3 + tmp31 + tmp50 + tmp51;  
                             EM_S[INDEX2(1,1,4)]+=tmp16 + tmp17 + tmp21 + tmp25 + tmp28 + tmp54 + tmp55;  
                             EM_S[INDEX2(1,2,4)]+=tmp10 + tmp6 + tmp7 + tmp8 + tmp9;  
                             EM_S[INDEX2(1,3,4)]+=tmp1 + tmp30 + tmp33 + tmp34 + tmp44 + tmp45;  
                             EM_S[INDEX2(2,0,4)]+=tmp1 + tmp34 + tmp38 + tmp40 + tmp42 + tmp43;  
                             EM_S[INDEX2(2,1,4)]+=tmp36 + tmp37 + tmp7 + tmp8 + tmp9;  
                             EM_S[INDEX2(2,2,4)]+=tmp15 + tmp19 + tmp22 + tmp23 + tmp28 + tmp48 + tmp49;  
                             EM_S[INDEX2(2,3,4)]+=tmp0 + tmp11 + tmp12 + tmp31 + tmp46 + tmp47;  
                             EM_S[INDEX2(3,0,4)]+=tmp27 + tmp28 + tmp29 + tmp7 + tmp8;  
                             EM_S[INDEX2(3,1,4)]+=tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35;  
                             EM_S[INDEX2(3,2,4)]+=tmp0 + tmp1 + tmp11 + tmp12 + tmp13 + tmp14;  
                             EM_S[INDEX2(3,3,4)]+=tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp9;  
                         } else { // constant data  
                             const double A_00 = A_p[INDEX2(0,0,2)];  
                             const double A_01 = A_p[INDEX2(0,1,2)];  
                             const double A_10 = A_p[INDEX2(1,0,2)];  
                             const double A_11 = A_p[INDEX2(1,1,2)];  
                             const double tmp0 = w9*(-A_01 - A_10);  
                             const double tmp1 = w9*(-A_01 + A_10);  
                             const double tmp2 = w9*(A_01 + A_10);  
                             const double tmp3 = w9*(A_01 - A_10);  
                             EM_S[INDEX2(0,0,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;  
                             EM_S[INDEX2(0,1,4)]+=8*A_00*w6 + 6*A_01*w1 + A_10*w9 + 4*A_11*w3;  
                             EM_S[INDEX2(0,2,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;  
                             EM_S[INDEX2(0,3,4)]+=4*A_00*w6 + A_11*w10 + tmp2;  
                             EM_S[INDEX2(1,0,4)]+=8*A_00*w6 - A_11*w10 + tmp3;  
                             EM_S[INDEX2(1,1,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;  
                             EM_S[INDEX2(1,2,4)]+=4*A_00*w6 + A_11*w10 + tmp0;  
                             EM_S[INDEX2(1,3,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;  
                             EM_S[INDEX2(2,0,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;  
                             EM_S[INDEX2(2,1,4)]+=4*A_00*w6 + A_11*w10 + tmp0;  
                             EM_S[INDEX2(2,2,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;  
                             EM_S[INDEX2(2,3,4)]+=8*A_00*w6 - A_11*w10 + tmp3;  
                             EM_S[INDEX2(3,0,4)]+=4*A_00*w6 + A_11*w10 + tmp2;  
                             EM_S[INDEX2(3,1,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;  
                             EM_S[INDEX2(3,2,4)]+=8*A_00*w6 - A_11*w10 + tmp1;  
                             EM_S[INDEX2(3,3,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;  
                         }  
                     }  
                     ///////////////  
                     // 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 = w15*(B_1_2 + B_1_3);  
                             const double tmp1 = w12*(B_1_0 + B_1_1);  
                             const double tmp2 = w15*(B_1_0 + B_1_1);  
                             const double tmp3 = w16*(-B_0_1 - B_0_3);  
                             const double tmp4 = w11*(-B_0_0 - B_0_2);  
                             const double tmp5 = w12*(B_1_2 + B_1_3);  
                             const double tmp6 = w15*(-B_1_0 - B_1_1);  
                             const double tmp7 = w12*(-B_1_2 - B_1_3);  
                             const double tmp8 = w15*(-B_1_2 - B_1_3);  
                             const double tmp9 = w12*(-B_1_0 - B_1_1);  
                             const double tmp10 = w11*(-B_0_1 - B_0_3);  
                             const double tmp11 = w16*(-B_0_0 - B_0_2);  
                             const double tmp12 = w16*(B_0_0 + B_0_2);  
                             const double tmp13 = w11*(B_0_1 + B_0_3);  
                             const double tmp14 = w11*(B_0_0 + B_0_2);  
                             const double tmp15 = w16*(B_0_1 + B_0_3);  
                             EM_S[INDEX2(0,0,4)]+=B_0_0*w13 + B_0_1*w11 + B_0_2*w16 + B_0_3*w14 + B_1_0*w17 + B_1_1*w15 + B_1_2*w12 + B_1_3*w18;  
                             EM_S[INDEX2(0,1,4)]+=B_0_0*w11 + B_0_1*w13 + B_0_2*w14 + B_0_3*w16 + tmp0 + tmp1;  
                             EM_S[INDEX2(0,2,4)]+=B_1_0*w12 + B_1_1*w18 + B_1_2*w17 + B_1_3*w15 + tmp14 + tmp15;  
                             EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp2 + tmp5;  
                             EM_S[INDEX2(1,0,4)]+=-B_0_0*w13 - B_0_1*w11 - B_0_2*w16 - B_0_3*w14 + tmp0 + tmp1;  
                             EM_S[INDEX2(1,1,4)]+=-B_0_0*w11 - B_0_1*w13 - B_0_2*w14 - B_0_3*w16 + B_1_0*w15 + B_1_1*w17 + B_1_2*w18 + B_1_3*w12;  
                             EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;  
                             EM_S[INDEX2(1,3,4)]+=B_1_0*w18 + B_1_1*w12 + B_1_2*w15 + B_1_3*w17 + tmp10 + tmp11;  
                             EM_S[INDEX2(2,0,4)]+=-B_1_0*w17 - B_1_1*w15 - B_1_2*w12 - B_1_3*w18 + tmp14 + tmp15;  
                             EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;  
                             EM_S[INDEX2(2,2,4)]+=B_0_0*w16 + B_0_1*w14 + B_0_2*w13 + B_0_3*w11 - B_1_0*w12 - B_1_1*w18 - B_1_2*w17 - B_1_3*w15;  
                             EM_S[INDEX2(2,3,4)]+=B_0_0*w14 + B_0_1*w16 + B_0_2*w11 + B_0_3*w13 + tmp6 + tmp7;  
                             EM_S[INDEX2(3,0,4)]+=tmp3 + tmp4 + tmp8 + tmp9;  
                             EM_S[INDEX2(3,1,4)]+=-B_1_0*w15 - B_1_1*w17 - B_1_2*w18 - B_1_3*w12 + tmp10 + tmp11;  
                             EM_S[INDEX2(3,2,4)]+=-B_0_0*w16 - B_0_1*w14 - B_0_2*w13 - B_0_3*w11 + tmp6 + tmp7;  
                             EM_S[INDEX2(3,3,4)]+=-B_0_0*w14 - B_0_1*w16 - B_0_2*w11 - B_0_3*w13 - B_1_0*w18 - B_1_1*w12 - B_1_2*w15 - B_1_3*w17;  
                         } else { // constant data  
                             const double B_0 = B_p[0];  
                             const double B_1 = B_p[1];  
                             EM_S[INDEX2(0,0,4)]+=B_0*w19 + 2*B_1*w20;  
                             EM_S[INDEX2(0,1,4)]+=B_0*w19 + B_1*w20;  
                             EM_S[INDEX2(0,2,4)]+=B_0*w19/2 + 2*B_1*w20;  
                             EM_S[INDEX2(0,3,4)]+=B_0*w19/2 + B_1*w20;  
                             EM_S[INDEX2(1,0,4)]+=-B_0*w19 + B_1*w20;  
                             EM_S[INDEX2(1,1,4)]+=-B_0*w19 + 2*B_1*w20;  
                             EM_S[INDEX2(1,2,4)]+=-B_0*w19/2 + B_1*w20;  
                             EM_S[INDEX2(1,3,4)]+=-B_0*w19/2 + 2*B_1*w20;  
                             EM_S[INDEX2(2,0,4)]+=B_0*w19/2 - 2*B_1*w20;  
                             EM_S[INDEX2(2,1,4)]+=B_0*w19/2 - B_1*w20;  
                             EM_S[INDEX2(2,2,4)]+=B_0*w19 - 2*B_1*w20;  
                             EM_S[INDEX2(2,3,4)]+=B_0*w19 - B_1*w20;  
                             EM_S[INDEX2(3,0,4)]+=-B_0*w19/2 - B_1*w20;  
                             EM_S[INDEX2(3,1,4)]+=-B_0*w19/2 - 2*B_1*w20;  
                             EM_S[INDEX2(3,2,4)]+=-B_0*w19 - B_1*w20;  
                             EM_S[INDEX2(3,3,4)]+=-B_0*w19 - 2*B_1*w20;  
                         }  
                     }  
                     ///////////////  
                     // 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 = w15*(C_1_2 + C_1_3);  
                             const double tmp1 = w12*(C_1_0 + C_1_1);  
                             const double tmp2 = w15*(-C_1_2 - C_1_3);  
                             const double tmp3 = w16*(C_0_0 + C_0_2);  
                             const double tmp4 = w11*(C_0_1 + C_0_3);  
                             const double tmp5 = w12*(-C_1_0 - C_1_1);  
                             const double tmp6 = w15*(-C_1_0 - C_1_1);  
                             const double tmp7 = w12*(-C_1_2 - C_1_3);  
                             const double tmp8 = w15*(C_1_0 + C_1_1);  
                             const double tmp9 = w12*(C_1_2 + C_1_3);  
                             const double tmp10 = w11*(-C_0_1 - C_0_3);  
                             const double tmp11 = w16*(-C_0_0 - C_0_2);  
                             const double tmp12 = w16*(-C_0_1 - C_0_3);  
                             const double tmp13 = w11*(-C_0_0 - C_0_2);  
                             const double tmp14 = w11*(C_0_0 + C_0_2);  
                             const double tmp15 = w16*(C_0_1 + C_0_3);  
                             EM_S[INDEX2(0,0,4)]+=C_0_0*w13 + C_0_1*w11 + C_0_2*w16 + C_0_3*w14 + C_1_0*w17 + C_1_1*w15 + C_1_2*w12 + C_1_3*w18;  
                             EM_S[INDEX2(0,1,4)]+=-C_0_0*w13 - C_0_1*w11 - C_0_2*w16 - C_0_3*w14 + tmp0 + tmp1;  
                             EM_S[INDEX2(0,2,4)]+=-C_1_0*w17 - C_1_1*w15 - C_1_2*w12 - C_1_3*w18 + tmp14 + tmp15;  
                             EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp2 + tmp5;  
                             EM_S[INDEX2(1,0,4)]+=C_0_0*w11 + C_0_1*w13 + C_0_2*w14 + C_0_3*w16 + tmp0 + tmp1;  
                             EM_S[INDEX2(1,1,4)]+=-C_0_0*w11 - C_0_1*w13 - C_0_2*w14 - C_0_3*w16 + C_1_0*w15 + C_1_1*w17 + C_1_2*w18 + C_1_3*w12;  
                             EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;  
                             EM_S[INDEX2(1,3,4)]+=-C_1_0*w15 - C_1_1*w17 - C_1_2*w18 - C_1_3*w12 + tmp10 + tmp11;  
                             EM_S[INDEX2(2,0,4)]+=C_1_0*w12 + C_1_1*w18 + C_1_2*w17 + C_1_3*w15 + tmp14 + tmp15;  
                             EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;  
                             EM_S[INDEX2(2,2,4)]+=C_0_0*w16 + C_0_1*w14 + C_0_2*w13 + C_0_3*w11 - C_1_0*w12 - C_1_1*w18 - C_1_2*w17 - C_1_3*w15;  
                             EM_S[INDEX2(2,3,4)]+=-C_0_0*w16 - C_0_1*w14 - C_0_2*w13 - C_0_3*w11 + tmp6 + tmp7;  
                             EM_S[INDEX2(3,0,4)]+=tmp3 + tmp4 + tmp8 + tmp9;  
                             EM_S[INDEX2(3,1,4)]+=C_1_0*w18 + C_1_1*w12 + C_1_2*w15 + C_1_3*w17 + tmp10 + tmp11;  
                             EM_S[INDEX2(3,2,4)]+=C_0_0*w14 + C_0_1*w16 + C_0_2*w11 + C_0_3*w13 + tmp6 + tmp7;  
                             EM_S[INDEX2(3,3,4)]+=-C_0_0*w14 - C_0_1*w16 - C_0_2*w11 - C_0_3*w13 - C_1_0*w18 - C_1_1*w12 - C_1_2*w15 - C_1_3*w17;  
                         } else { // constant data  
                             const double C_0 = C_p[0];  
                             const double C_1 = C_p[1];  
                             EM_S[INDEX2(0,0,4)]+=C_0*w19 + 2*C_1*w20;  
                             EM_S[INDEX2(0,1,4)]+=-C_0*w19 + C_1*w20;  
                             EM_S[INDEX2(0,2,4)]+=C_0*w19/2 - 2*C_1*w20;  
                             EM_S[INDEX2(0,3,4)]+=-C_0*w19/2 - C_1*w20;  
                             EM_S[INDEX2(1,0,4)]+=C_0*w19 + C_1*w20;  
                             EM_S[INDEX2(1,1,4)]+=-C_0*w19 + 2*C_1*w20;  
                             EM_S[INDEX2(1,2,4)]+=C_0*w19/2 - C_1*w20;  
                             EM_S[INDEX2(1,3,4)]+=-C_0*w19/2 - 2*C_1*w20;  
                             EM_S[INDEX2(2,0,4)]+=C_0*w19/2 + 2*C_1*w20;  
                             EM_S[INDEX2(2,1,4)]+=-C_0*w19/2 + C_1*w20;  
                             EM_S[INDEX2(2,2,4)]+=C_0*w19 - 2*C_1*w20;  
                             EM_S[INDEX2(2,3,4)]+=-C_0*w19 - C_1*w20;  
                             EM_S[INDEX2(3,0,4)]+=C_0*w19/2 + C_1*w20;  
                             EM_S[INDEX2(3,1,4)]+=-C_0*w19/2 + 2*C_1*w20;  
                             EM_S[INDEX2(3,2,4)]+=C_0*w19 - C_1*w20;  
                             EM_S[INDEX2(3,3,4)]+=-C_0*w19 - 2*C_1*w20;  
                         }  
                     }  
                     ///////////////  
                     // 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 = w21*(D_0 + D_1);  
                             const double tmp1 = w22*(D_2 + D_3);  
                             const double tmp2 = w23*(D_0 + D_1 + D_2 + D_3);  
                             const double tmp3 = w21*(D_2 + D_3);  
                             const double tmp4 = w22*(D_0 + D_1);  
                             const double tmp5 = w23*(D_1 + D_2);  
                             const double tmp6 = w21*(D_1 + D_3);  
                             const double tmp7 = w22*(D_0 + D_2);  
                             const double tmp8 = w21*(D_0 + D_2);  
                             const double tmp9 = w22*(D_1 + D_3);  
                             const double tmp10 = w23*(D_0 + D_3);  
                             EM_S[INDEX2(0,0,4)]+=D_0*w24 + D_3*w25 + tmp5;  
                             EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1;  
                             EM_S[INDEX2(0,2,4)]+=tmp8 + tmp9;  
                             EM_S[INDEX2(0,3,4)]+=tmp2;  
                             EM_S[INDEX2(1,0,4)]+=tmp0 + tmp1;  
                             EM_S[INDEX2(1,1,4)]+=D_1*w24 + D_2*w25 + tmp10;  
                             EM_S[INDEX2(1,2,4)]+=tmp2;  
                             EM_S[INDEX2(1,3,4)]+=tmp6 + tmp7;  
                             EM_S[INDEX2(2,0,4)]+=tmp8 + tmp9;  
                             EM_S[INDEX2(2,1,4)]+=tmp2;  
                             EM_S[INDEX2(2,2,4)]+=D_1*w25 + D_2*w24 + tmp10;  
                             EM_S[INDEX2(2,3,4)]+=tmp3 + tmp4;  
                             EM_S[INDEX2(3,0,4)]+=tmp2;  
                             EM_S[INDEX2(3,1,4)]+=tmp6 + tmp7;  
                             EM_S[INDEX2(3,2,4)]+=tmp3 + tmp4;  
                             EM_S[INDEX2(3,3,4)]+=D_0*w25 + D_3*w24 + tmp5;  
                         } else { // constant data  
                             const double D_0 = D_p[0];  
                             EM_S[INDEX2(0,0,4)]+=16*D_0*w23;  
                             EM_S[INDEX2(0,1,4)]+=8*D_0*w23;  
                             EM_S[INDEX2(0,2,4)]+=8*D_0*w23;  
                             EM_S[INDEX2(0,3,4)]+=4*D_0*w23;  
                             EM_S[INDEX2(1,0,4)]+=8*D_0*w23;  
                             EM_S[INDEX2(1,1,4)]+=16*D_0*w23;  
                             EM_S[INDEX2(1,2,4)]+=4*D_0*w23;  
                             EM_S[INDEX2(1,3,4)]+=8*D_0*w23;  
                             EM_S[INDEX2(2,0,4)]+=8*D_0*w23;  
                             EM_S[INDEX2(2,1,4)]+=4*D_0*w23;  
                             EM_S[INDEX2(2,2,4)]+=16*D_0*w23;  
                             EM_S[INDEX2(2,3,4)]+=8*D_0*w23;  
                             EM_S[INDEX2(3,0,4)]+=4*D_0*w23;  
                             EM_S[INDEX2(3,1,4)]+=8*D_0*w23;  
                             EM_S[INDEX2(3,2,4)]+=8*D_0*w23;  
                             EM_S[INDEX2(3,3,4)]+=16*D_0*w23;  
                         }  
                     }  
                     ///////////////  
                     // 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 = 6*w15*(X_1_1 + X_1_3);  
                             const double tmp1 = 6*w16*(X_0_2 + X_0_3);  
                             const double tmp2 = 6*w11*(X_0_0 + X_0_1);  
                             const double tmp3 = 6*w12*(X_1_0 + X_1_2);  
                             const double tmp4 = 6*w15*(X_1_0 + X_1_2);  
                             const double tmp5 = w27*(X_0_2 + X_0_3);  
                             const double tmp6 = w26*(X_0_0 + X_0_1);  
                             const double tmp7 = 6*w12*(X_1_1 + X_1_3);  
                             const double tmp8 = w29*(X_1_1 + X_1_3);  
                             const double tmp9 = w27*(-X_0_0 - X_0_1);  
                             const double tmp10 = w28*(X_1_0 + X_1_2);  
                             const double tmp11 = w26*(-X_0_2 - X_0_3);  
                             const double tmp12 = w29*(X_1_0 + X_1_2);  
                             const double tmp13 = w27*(X_0_0 + X_0_1);  
                             const double tmp14 = w28*(X_1_1 + X_1_3);  
                             const double tmp15 = w26*(X_0_2 + X_0_3);  
                             EM_F[0]+=tmp0 + tmp1 + tmp2 + tmp3;  
                             EM_F[1]+=tmp4 + tmp5 + tmp6 + tmp7;  
                             EM_F[2]+=tmp10 + tmp11 + tmp8 + tmp9;  
                             EM_F[3]+=tmp12 + tmp13 + tmp14 + tmp15;  
                         } else { // constant data  
                             const double X_0 = X_p[0];  
                             const double X_1 = X_p[1];  
                             EM_F[0]+=3*X_0*w19 + 6*X_1*w20;  
                             EM_F[1]+=-3*X_0*w19 + 6*X_1*w20;  
                             EM_F[2]+=3*X_0*w19 - 6*X_1*w20;  
                             EM_F[3]+=-3*X_0*w19 - 6*X_1*w20;  
                         }  
                     }  
                     ///////////////  
                     // 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 = 6*w23*(Y_1 + Y_2);  
                             const double tmp1 = 6*w23*(Y_0 + Y_3);  
                             EM_F[0]+=6*Y_0*w21 + 6*Y_3*w22 + tmp0;  
                             EM_F[1]+=6*Y_1*w21 + 6*Y_2*w22 + tmp1;  
                             EM_F[2]+=6*Y_1*w22 + 6*Y_2*w21 + tmp1;  
                             EM_F[3]+=6*Y_0*w22 + 6*Y_3*w21 + tmp0;  
                         } else { // constant data  
                             const double Y_0 = Y_p[0];  
                             EM_F[0]+=36*Y_0*w23;  
                             EM_F[1]+=36*Y_0*w23;  
                             EM_F[2]+=36*Y_0*w23;  
                             EM_F[3]+=36*Y_0*w23;  
                         }  
                     }  
                     /* GENERATOR SNIP_PDE_SINGLE BOTTOM */  
1870    
1871                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)  namespace
1872                      const index_t firstNode=m_NN[0]*k1+k0;  {
1873                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);      // Calculates a guassian blur colvolution matrix for 2D
1874                  } // end k0 loop      // See wiki article on the subject    
1875              } // end k1 loop      double* get2DGauss(unsigned radius, double sigma)
1876          } // end of colouring      {
1877      } // end of parallel region          double* arr=new double[(radius*2+1)*(radius*2+1)];
1878            double common=M_1_PI*0.5*1/(sigma*sigma);
1879            double total=0;
1880            int r=static_cast<int>(radius);
1881            for (int y=-r;y<=r;++y)
1882            {
1883                for (int x=-r;x<=r;++x)
1884                {        
1885                    arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));
1886                    total+=arr[(x+r)+(y+r)*(r*2+1)];
1887                }
1888            }
1889            double invtotal=1/total;
1890            for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)
1891            {
1892                arr[p]*=invtotal;
1893            }
1894            return arr;
1895        }
1896        
1897        // applies conv to source to get a point.
1898        // (xp, yp) are the coords in the source matrix not the destination matrix
1899        double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)
1900        {
1901            size_t bx=xp-radius, by=yp-radius;
1902            size_t sbase=bx+by*width;
1903            double result=0;
1904            for (int y=0;y<2*radius+1;++y)
1905            {        
1906                for (int x=0;x<2*radius+1;++x)
1907                {
1908                    result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];
1909                }
1910            }
1911            return result;      
1912        }
1913  }  }
1914    
 //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 w0 = .25;  
     const double w1 = .125*m_dx[0];  
     const double w2 = .125*m_dx[1];  
     const double w3 = .0625*m_dx[0]*m_dx[1];  
     const double w4 = .25*m_dx[0]/m_dx[1];  
     const double w5 = .25*m_dx[1]/m_dx[0];  
1915    
1916      rhs.requireWrite();  /* This is a wrapper for filtered (and non-filtered) randoms
1917  #pragma omp parallel   * For detailed doco see randomFillWorker
1918    */
1919    escript::Data Rectangle::randomFill(const escript::DataTypes::ShapeType& shape,
1920           const escript::FunctionSpace& what,
1921           long seed, const boost::python::tuple& filter) const
1922    {
1923        int numvals=escript::DataTypes::noValues(shape);
1924        if (len(filter)>0 && (numvals!=1))
1925      {      {
1926          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          throw RipleyException("Ripley only supports filters for scalar data.");
1927  #pragma omp for      }
1928              for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {      escript::Data res=randomFillWorker(shape, seed, filter);
1929                  for (index_t k0=0; k0<m_NE[0]; ++k0)  {      if (res.getFunctionSpace()!=what)
1930                      bool add_EM_S=false;      {
1931                      bool add_EM_F=false;          escript::Data r=escript::Data(res, what);
1932                      vector<double> EM_S(4*4, 0);          return r;
1933                      vector<double> EM_F(4, 0);      }
1934                      const index_t e = k0 + m_NE[0]*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 = (A_01 + A_10)*w0;  
                         const double tmp1 = A_00*w5;  
                         const double tmp2 = A_01*w0;  
                         const double tmp3 = A_10*w0;  
                         const double tmp4 = A_11*w4;  
                         EM_S[INDEX2(0,0,4)]+=tmp4 + tmp0 + tmp1;  
                         EM_S[INDEX2(1,0,4)]+=tmp4 - tmp1 + tmp3 - tmp2;  
                         EM_S[INDEX2(2,0,4)]+=tmp2 - tmp3 - tmp4 + tmp1;  
                         EM_S[INDEX2(3,0,4)]+=-tmp1 - tmp4 - tmp0;  
                         EM_S[INDEX2(0,1,4)]+=tmp4 - tmp1 + tmp2 - tmp3;  
                         EM_S[INDEX2(1,1,4)]+=tmp4 + tmp1 - tmp0;  
                         EM_S[INDEX2(2,1,4)]+=-tmp1 + tmp0 - tmp4;  
                         EM_S[INDEX2(3,1,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;  
                         EM_S[INDEX2(0,2,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;  
                         EM_S[INDEX2(1,2,4)]+=-tmp1 + tmp0 - tmp4;  
                         EM_S[INDEX2(2,2,4)]+=tmp4 + tmp1 - tmp0;  
                         EM_S[INDEX2(3,2,4)]+=tmp4 - tmp1 + tmp2 - tmp3;  
                         EM_S[INDEX2(0,3,4)]+=-tmp1 - tmp4 - tmp0;  
                         EM_S[INDEX2(1,3,4)]+=tmp2 - tmp3 - tmp4 + tmp1;  
                         EM_S[INDEX2(2,3,4)]+=tmp4 - tmp1 + tmp3 - tmp2;  
                         EM_S[INDEX2(3,3,4)]+=tmp4 + tmp0 + tmp1;  
                     }  
                     ///////////////  
                     // process B //  
                     ///////////////  
                     if (!B.isEmpty()) {  
                         add_EM_S=true;  
                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);  
                         const double tmp0 = B_p[0]*w2;  
                         const double tmp1 = B_p[1]*w1;  
                         EM_S[INDEX2(0,0,4)]+=-tmp0 - tmp1;  
                         EM_S[INDEX2(1,0,4)]+= tmp0 - tmp1;  
                         EM_S[INDEX2(2,0,4)]+= tmp1 - tmp0;  
                         EM_S[INDEX2(3,0,4)]+= tmp0 + tmp1;  
                         EM_S[INDEX2(0,1,4)]+=-tmp0 - tmp1;  
                         EM_S[INDEX2(1,1,4)]+= tmp0 - tmp1;  
                         EM_S[INDEX2(2,1,4)]+= tmp1 - tmp0;  
                         EM_S[INDEX2(3,1,4)]+= tmp0 + tmp1;  
                         EM_S[INDEX2(0,2,4)]+=-tmp0 - tmp1;  
                         EM_S[INDEX2(1,2,4)]+= tmp0 - tmp1;  
                         EM_S[INDEX2(2,2,4)]+= tmp1 - tmp0;  
                         EM_S[INDEX2(3,2,4)]+= tmp0 + tmp1;  
                         EM_S[INDEX2(0,3,4)]+=-tmp0 - tmp1;  
                         EM_S[INDEX2(1,3,4)]+= tmp0 - tmp1;  
                         EM_S[INDEX2(2,3,4)]+= tmp1 - tmp0;  
                         EM_S[INDEX2(3,3,4)]+= tmp0 + tmp1;  
                     }  
                     ///////////////  
                     // process C //  
                     ///////////////  
                     if (!C.isEmpty()) {  
                         add_EM_S=true;  
                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);  
                         const double tmp0 = C_p[0]*w2;  
                         const double tmp1 = C_p[1]*w1;  
                         EM_S[INDEX2(0,0,4)]+=-tmp1 - tmp0;  
                         EM_S[INDEX2(1,0,4)]+=-tmp1 - tmp0;  
                         EM_S[INDEX2(2,0,4)]+=-tmp1 - tmp0;  
                         EM_S[INDEX2(3,0,4)]+=-tmp1 - tmp0;  
                         EM_S[INDEX2(0,1,4)]+= tmp0 - tmp1;  
                         EM_S[INDEX2(1,1,4)]+= tmp0 - tmp1;  
                         EM_S[INDEX2(2,1,4)]+= tmp0 - tmp1;  
                         EM_S[INDEX2(3,1,4)]+= tmp0 - tmp1;  
                         EM_S[INDEX2(0,2,4)]+= tmp1 - tmp0;  
                         EM_S[INDEX2(1,2,4)]+= tmp1 - tmp0;  
                         EM_S[INDEX2(2,2,4)]+= tmp1 - tmp0;  
                         EM_S[INDEX2(3,2,4)]+= tmp1 - tmp0;  
                         EM_S[INDEX2(0,3,4)]+= tmp0 + tmp1;  
                         EM_S[INDEX2(1,3,4)]+= tmp0 + tmp1;  
                         EM_S[INDEX2(2,3,4)]+= tmp0 + tmp1;  
                         EM_S[INDEX2(3,3,4)]+= tmp0 + tmp1;  
                     }  
                     ///////////////  
                     // process D //  
                     ///////////////  
                     if (!D.isEmpty()) {  
                         add_EM_S=true;  
                         const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);  
                         EM_S[INDEX2(0,0,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(1,0,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(2,0,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(3,0,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(0,1,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(1,1,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(2,1,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(3,1,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(0,2,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(1,2,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(2,2,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(3,2,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(0,3,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(1,3,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(2,3,4)]+=D_p[0]*w3;  
                         EM_S[INDEX2(3,3,4)]+=D_p[0]*w3;  
                     }  
                     ///////////////  
                     // process X //  
                     ///////////////  
                     if (!X.isEmpty()) {  
                         add_EM_F=true;  
                         const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);  
                         const double tmp0 = 4*X_p[0]*w2;  
                         const double tmp1 = 4*X_p[1]*w1;  
                         EM_F[0]+=-tmp0 - tmp1;  
                         EM_F[1]+=-tmp1 + tmp0;  
                         EM_F[2]+=-tmp0 + tmp1;  
                         EM_F[3]+= tmp0 + tmp1;  
                     }  
                     ///////////////  
                     // process Y //  
                     ///////////////  
                     if (!Y.isEmpty()) {  
                         add_EM_F=true;  
                         const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);  
                         EM_F[0]+=4*Y_p[0]*w3;  
                         EM_F[1]+=4*Y_p[0]*w3;  
                         EM_F[2]+=4*Y_p[0]*w3;  
                         EM_F[3]+=4*Y_p[0]*w3;  
                     }  
   
                     // add to matrix (if add_EM_S) and RHS (if add_EM_F)  
                     const index_t firstNode=m_NN[0]*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  
1935  }  }
1936    
 //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  
 {  
     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;  
     }  
   
     /* GENERATOR SNIP_PDE_SYSTEM_PRE TOP */  
     const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];  
     const double w1 = 0.041666666666666666667;  
     const double w2 = -0.15550211698203655390;  
     const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];  
     const double w4 = -0.01116454968463011277*m_dx[1]/m_dx[0];  
     const double w5 = 0.011164549684630112770;  
     const double w6 = -0.041666666666666666667*m_dx[1]/m_dx[0];  
     const double w7 = 0.1555021169820365539*m_dx[0]/m_dx[1];  
     const double w8 = 0.01116454968463011277*m_dx[0]/m_dx[1];  
     const double w9 = -0.25;  
     const double w10 = -0.16666666666666666667*m_dx[0]/m_dx[1];  
     const double w11 = -0.032861463941450536761*m_dx[1];  
     const double w12 = -0.032861463941450536761*m_dx[0];  
     const double w13 = -0.12264065304058601714*m_dx[1];  
     const double w14 = -0.0023593469594139828636*m_dx[1];  
     const double w15 = -0.008805202725216129906*m_dx[0];  
     const double w16 = -0.008805202725216129906*m_dx[1];  
     const double w17 = -0.12264065304058601714*m_dx[0];  
     const double w18 = -0.0023593469594139828636*m_dx[0];  
     const double w19 = -0.16666666666666666667*m_dx[1];  
     const double w20 = -0.083333333333333333333*m_dx[0];  
     const double w21 = 0.025917019497006092316*m_dx[0]*m_dx[1];  
     const double w22 = 0.0018607582807716854616*m_dx[0]*m_dx[1];  
     const double w23 = 0.0069444444444444444444*m_dx[0]*m_dx[1];  
     const double w24 = 0.09672363354357992482*m_dx[0]*m_dx[1];  
     const double w25 = 0.000498588678642297402*m_dx[0]*m_dx[1];  
     const double w26 = 0.19716878364870322056*m_dx[1];  
     const double w27 = 0.052831216351296779436*m_dx[1];  
     const double w28 = 0.19716878364870322056*m_dx[0];  
     const double w29 = 0.052831216351296779436*m_dx[0];  
     /* GENERATOR SNIP_PDE_SYSTEM_PRE BOTTOM */  
1937    
1938      rhs.requireWrite();  /* This routine produces a Data object filled with smoothed random data.
1939  #pragma omp parallel  The dimensions of the rectangle being filled are internal[0] x internal[1] points.
1940      {  A parameter radius  gives the size of the stencil used for the smoothing.
1941          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
1942  #pragma omp for  in order to complete the stencil.
             for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {  
                 for (index_t k0=0; k0<m_NE[0]; ++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_NE[0]*k1;  
                     /* GENERATOR SNIP_PDE_SYSTEM TOP */  
                     ///////////////  
                     // 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 = w3*(A_11_0 + A_11_1 + A_11_2 + A_11_3);  
                                     const double tmp1 = w1*(A_01_0 + A_01_3 - A_10_1 - A_10_2);  
                                     const double tmp2 = w4*(A_00_2 + A_00_3);  
                                     const double tmp3 = w0*(A_00_0 + A_00_1);  
                                     const double tmp4 = w5*(A_01_2 - A_10_3);  
                                     const double tmp5 = w2*(-A_01_1 + A_10_0);  
                                     const double tmp6 = w5*(A_01_3 + A_10_0);  
                                     const double tmp7 = w3*(-A_11_0 - A_11_1 - A_11_2 - A_11_3);  
                                     const double tmp8 = w6*(A_00_0 + A_00_1 + A_00_2 + A_00_3);  
                                     const double tmp9 = w1*(A_01_1 + A_01_2 + A_10_1 + A_10_2);  
                                     const double tmp10 = w2*(-A_01_0 - A_10_3);  
                                     const double tmp11 = w4*(A_00_0 + A_00_1);  
                                     const double tmp12 = w0*(A_00_2 + A_00_3);  
                                     const double tmp13 = w5*(A_01_1 - A_10_0);  
                                     const double tmp14 = w2*(-A_01_2 + A_10_3);  
                                     const double tmp15 = w7*(A_11_0 + A_11_2);  
                                     const double tmp16 = w4*(-A_00_2 - A_00_3);  
                                     const double tmp17 = w0*(-A_00_0 - A_00_1);  
                                     const double tmp18 = w5*(A_01_3 + A_10_3);  
                                     const double tmp19 = w8*(A_11_1 + A_11_3);  
                                     const double tmp20 = w2*(-A_01_0 - A_10_0);  
                                     const double tmp21 = w7*(A_11_1 + A_11_3);  
                                     const double tmp22 = w4*(-A_00_0 - A_00_1);  
                                     const double tmp23 = w0*(-A_00_2 - A_00_3);  
                                     const double tmp24 = w5*(A_01_0 + A_10_0);  
                                     const double tmp25 = w8*(A_11_0 + A_11_2);  
                                     const double tmp26 = w2*(-A_01_3 - A_10_3);  
                                     const double tmp27 = w5*(-A_01_1 - A_10_2);  
                                     const double tmp28 = w1*(-A_01_0 - A_01_3 - A_10_0 - A_10_3);  
                                     const double tmp29 = w2*(A_01_2 + A_10_1);  
                                     const double tmp30 = w7*(-A_11_1 - A_11_3);  
                                     const double tmp31 = w1*(-A_01_1 - A_01_2 + A_10_0 + A_10_3);  
                                     const double tmp32 = w5*(-A_01_0 + A_10_2);  
                                     const double tmp33 = w8*(-A_11_0 - A_11_2);  
                                     const double tmp34 = w6*(-A_00_0 - A_00_1 - A_00_2 - A_00_3);  
                                     const double tmp35 = w2*(A_01_3 - A_10_1);  
                                     const double tmp36 = w5*(A_01_0 + A_10_3);  
                                     const double tmp37 = w2*(-A_01_3 - A_10_0);  
                                     const double tmp38 = w7*(-A_11_0 - A_11_2);  
                                     const double tmp39 = w5*(-A_01_3 + A_10_1);  
                                     const double tmp40 = w8*(-A_11_1 - A_11_3);  
                                     const double tmp41 = w2*(A_01_0 - A_10_2);  
                                     const double tmp42 = w5*(A_01_1 - A_10_3);  
                                     const double tmp43 = w2*(-A_01_2 + A_10_0);  
                                     const double tmp44 = w5*(A_01_2 - A_10_0);  
                                     const double tmp45 = w2*(-A_01_1 + A_10_3);  
                                     const double tmp46 = w5*(-A_01_0 + A_10_1);  
                                     const double tmp47 = w2*(A_01_3 - A_10_2);  
                                     const double tmp48 = w5*(-A_01_1 - A_10_1);  
                                     const double tmp49 = w2*(A_01_2 + A_10_2);  
                                     const double tmp50 = w5*(-A_01_3 + A_10_2);  
                                     const double tmp51 = w2*(A_01_0 - A_10_1);  
                                     const double tmp52 = w5*(-A_01_2 - A_10_1);  
                                     const double tmp53 = w2*(A_01_1 + A_10_2);  
                                     const double tmp54 = w5*(-A_01_2 - A_10_2);  
                                     const double tmp55 = w2*(A_01_1 + A_10_1);  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp9;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp31 + tmp34 + tmp38 + tmp39 + tmp40 + tmp41;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp28 + tmp52 + tmp53 + tmp7 + tmp8;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp2 + tmp3 + tmp31 + tmp50 + tmp51;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp16 + tmp17 + tmp21 + tmp25 + tmp28 + tmp54 + tmp55;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10 + tmp6 + tmp7 + tmp8 + tmp9;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1 + tmp30 + tmp33 + tmp34 + tmp44 + tmp45;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1 + tmp34 + tmp38 + tmp40 + tmp42 + tmp43;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp36 + tmp37 + tmp7 + tmp8 + tmp9;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp15 + tmp19 + tmp22 + tmp23 + tmp28 + tmp48 + tmp49;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0 + tmp11 + tmp12 + tmp31 + tmp46 + tmp47;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp27 + tmp28 + tmp29 + tmp7 + tmp8;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0 + tmp1 + tmp11 + tmp12 + tmp13 + tmp14;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp9;  
                                 }  
                             }  
                         } 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 = w9*(-A_01 - A_10);  
                                     const double tmp1 = w9*(-A_01 + A_10);  
                                     const double tmp2 = w9*(A_01 + A_10);  
                                     const double tmp3 = w9*(A_01 - A_10);  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*A_00*w6 + 6*A_01*w1 + A_10*w9 + 4*A_11*w3;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp2;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp3;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp0;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp0;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp3;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp2;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // 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 = w15*(B_1_2 + B_1_3);  
                                     const double tmp1 = w12*(B_1_0 + B_1_1);  
                                     const double tmp2 = w15*(B_1_0 + B_1_1);  
                                     const double tmp3 = w16*(-B_0_1 - B_0_3);  
                                     const double tmp4 = w11*(-B_0_0 - B_0_2);  
                                     const double tmp5 = w12*(B_1_2 + B_1_3);  
                                     const double tmp6 = w15*(-B_1_0 - B_1_1);  
                                     const double tmp7 = w12*(-B_1_2 - B_1_3);  
                                     const double tmp8 = w15*(-B_1_2 - B_1_3);  
                                     const double tmp9 = w12*(-B_1_0 - B_1_1);  
                                     const double tmp10 = w11*(-B_0_1 - B_0_3);  
                                     const double tmp11 = w16*(-B_0_0 - B_0_2);  
                                     const double tmp12 = w16*(B_0_0 + B_0_2);  
                                     const double tmp13 = w11*(B_0_1 + B_0_3);  
                                     const double tmp14 = w11*(B_0_0 + B_0_2);  
                                     const double tmp15 = w16*(B_0_1 + B_0_3);  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=B_0_0*w13 + B_0_1*w11 + B_0_2*w16 + B_0_3*w14 + B_1_0*w17 + B_1_1*w15 + B_1_2*w12 + B_1_3*w18;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=B_0_0*w11 + B_0_1*w13 + B_0_2*w14 + B_0_3*w16 + tmp0 + tmp1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=B_1_0*w12 + B_1_1*w18 + B_1_2*w17 + B_1_3*w15 + tmp14 + tmp15;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp2 + tmp5;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-B_0_0*w13 - B_0_1*w11 - B_0_2*w16 - B_0_3*w14 + tmp0 + tmp1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-B_0_0*w11 - B_0_1*w13 - B_0_2*w14 - B_0_3*w16 + B_1_0*w15 + B_1_1*w17 + B_1_2*w18 + B_1_3*w12;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=B_1_0*w18 + B_1_1*w12 + B_1_2*w15 + B_1_3*w17 + tmp10 + tmp11;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-B_1_0*w17 - B_1_1*w15 - B_1_2*w12 - B_1_3*w18 + tmp14 + tmp15;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=B_0_0*w16 + B_0_1*w14 + B_0_2*w13 + B_0_3*w11 - B_1_0*w12 - B_1_1*w18 - B_1_2*w17 - B_1_3*w15;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=B_0_0*w14 + B_0_1*w16 + B_0_2*w11 + B_0_3*w13 + tmp6 + tmp7;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp3 + tmp4 + tmp8 + tmp9;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-B_1_0*w15 - B_1_1*w17 - B_1_2*w18 - B_1_3*w12 + tmp10 + tmp11;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-B_0_0*w16 - B_0_1*w14 - B_0_2*w13 - B_0_3*w11 + tmp6 + tmp7;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-B_0_0*w14 - B_0_1*w16 - B_0_2*w11 - B_0_3*w13 - B_1_0*w18 - B_1_1*w12 - B_1_2*w15 - B_1_3*w17;  
                                 }  
                             }  
                         } 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)];  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=B_0*w19 + 2*B_1*w20;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=B_0*w19 + B_1*w20;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=B_0*w19/2 + 2*B_1*w20;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=B_0*w19/2 + B_1*w20;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-B_0*w19 + B_1*w20;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-B_0*w19 + 2*B_1*w20;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-B_0*w19/2 + B_1*w20;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-B_0*w19/2 + 2*B_1*w20;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=B_0*w19/2 - 2*B_1*w20;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=B_0*w19/2 - B_1*w20;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=B_0*w19 - 2*B_1*w20;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=B_0*w19 - B_1*w20;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-B_0*w19/2 - B_1*w20;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-B_0*w19/2 - 2*B_1*w20;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-B_0*w19 - B_1*w20;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-B_0*w19 - 2*B_1*w20;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // 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 tmp0 = w15*(C_1_2 + C_1_3);  
                                     const double tmp1 = w12*(C_1_0 + C_1_1);  
                                     const double tmp2 = w15*(-C_1_2 - C_1_3);  
                                     const double tmp3 = w16*(C_0_0 + C_0_2);  
                                     const double tmp4 = w11*(C_0_1 + C_0_3);  
                                     const double tmp5 = w12*(-C_1_0 - C_1_1);  
                                     const double tmp6 = w15*(-C_1_0 - C_1_1);  
                                     const double tmp7 = w12*(-C_1_2 - C_1_3);  
                                     const double tmp8 = w15*(C_1_0 + C_1_1);  
                                     const double tmp9 = w12*(C_1_2 + C_1_3);  
                                     const double tmp10 = w11*(-C_0_1 - C_0_3);  
                                     const double tmp11 = w16*(-C_0_0 - C_0_2);  
                                     const double tmp12 = w16*(-C_0_1 - C_0_3);  
                                     const double tmp13 = w11*(-C_0_0 - C_0_2);  
                                     const double tmp14 = w11*(C_0_0 + C_0_2);  
                                     const double tmp15 = w16*(C_0_1 + C_0_3);  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=C_0_0*w13 + C_0_1*w11 + C_0_2*w16 + C_0_3*w14 + C_1_0*w17 + C_1_1*w15 + C_1_2*w12 + C_1_3*w18;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-C_0_0*w13 - C_0_1*w11 - C_0_2*w16 - C_0_3*w14 + tmp0 + tmp1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-C_1_0*w17 - C_1_1*w15 - C_1_2*w12 - C_1_3*w18 + tmp14 + tmp15;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp2 + tmp5;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=C_0_0*w11 + C_0_1*w13 + C_0_2*w14 + C_0_3*w16 + tmp0 + tmp1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-C_0_0*w11 - C_0_1*w13 - C_0_2*w14 - C_0_3*w16 + C_1_0*w15 + C_1_1*w17 + C_1_2*w18 + C_1_3*w12;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-C_1_0*w15 - C_1_1*w17 - C_1_2*w18 - C_1_3*w12 + tmp10 + tmp11;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=C_1_0*w12 + C_1_1*w18 + C_1_2*w17 + C_1_3*w15 + tmp14 + tmp15;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=C_0_0*w16 + C_0_1*w14 + C_0_2*w13 + C_0_3*w11 - C_1_0*w12 - C_1_1*w18 - C_1_2*w17 - C_1_3*w15;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-C_0_0*w16 - C_0_1*w14 - C_0_2*w13 - C_0_3*w11 + tmp6 + tmp7;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp3 + tmp4 + tmp8 + tmp9;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=C_1_0*w18 + C_1_1*w12 + C_1_2*w15 + C_1_3*w17 + tmp10 + tmp11;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=C_0_0*w14 + C_0_1*w16 + C_0_2*w11 + C_0_3*w13 + tmp6 + tmp7;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-C_0_0*w14 - C_0_1*w16 - C_0_2*w11 - C_0_3*w13 - C_1_0*w18 - C_1_1*w12 - C_1_2*w15 - C_1_3*w17;  
                                 }  
                             }  
                         } 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)];  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=C_0*w19 + 2*C_1*w20;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-C_0*w19 + C_1*w20;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=C_0*w19/2 - 2*C_1*w20;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-C_0*w19/2 - C_1*w20;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=C_0*w19 + C_1*w20;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-C_0*w19 + 2*C_1*w20;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=C_0*w19/2 - C_1*w20;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-C_0*w19/2 - 2*C_1*w20;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=C_0*w19/2 + 2*C_1*w20;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-C_0*w19/2 + C_1*w20;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=C_0*w19 - 2*C_1*w20;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-C_0*w19 - C_1*w20;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=C_0*w19/2 + C_1*w20;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-C_0*w19/2 + 2*C_1*w20;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=C_0*w19 - C_1*w20;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-C_0*w19 - 2*C_1*w20;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // 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 tmp0 = w21*(D_0 + D_1);  
                                     const double tmp1 = w22*(D_2 + D_3);  
                                     const double tmp2 = w23*(D_0 + D_1 + D_2 + D_3);  
                                     const double tmp3 = w21*(D_2 + D_3);  
                                     const double tmp4 = w22*(D_0 + D_1);  
                                     const double tmp5 = w23*(D_1 + D_2);  
                                     const double tmp6 = w21*(D_1 + D_3);  
                                     const double tmp7 = w22*(D_0 + D_2);  
                                     const double tmp8 = w21*(D_0 + D_2);  
                                     const double tmp9 = w22*(D_1 + D_3);  
                                     const double tmp10 = w23*(D_0 + D_3);  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=D_0*w24 + D_3*w25 + tmp5;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp8 + tmp9;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=D_1*w24 + D_2*w25 + tmp10;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp6 + tmp7;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp8 + tmp9;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp2;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=D_1*w25 + D_2*w24 + tmp10;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp3 + tmp4;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp6 + tmp7;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3 + tmp4;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=D_0*w25 + D_3*w24 + tmp5;  
                                 }  
                              }  
                         } 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)];  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=16*D_0*w23;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*D_0*w23;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=8*D_0*w23;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*D_0*w23;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*D_0*w23;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=16*D_0*w23;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*D_0*w23;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=8*D_0*w23;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=8*D_0*w23;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*D_0*w23;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=16*D_0*w23;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*D_0*w23;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*D_0*w23;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=8*D_0*w23;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*D_0*w23;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=16*D_0*w23;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // 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 tmp0 = 6*w15*(X_1_1 + X_1_3);  
                                 const double tmp1 = 6*w16*(X_0_2 + X_0_3);  
                                 const double tmp2 = 6*w11*(X_0_0 + X_0_1);  
                                 const double tmp3 = 6*w12*(X_1_0 + X_1_2);  
                                 const double tmp4 = 6*w15*(X_1_0 + X_1_2);  
                                 const double tmp5 = w27*(X_0_2 + X_0_3);  
                                 const double tmp6 = w26*(X_0_0 + X_0_1);  
                                 const double tmp7 = 6*w12*(X_1_1 + X_1_3);  
                                 const double tmp8 = w29*(X_1_1 + X_1_3);  
                                 const double tmp9 = w27*(-X_0_0 - X_0_1);  
                                 const double tmp10 = w28*(X_1_0 + X_1_2);  
                                 const double tmp11 = w26*(-X_0_2 - X_0_3);  
                                 const double tmp12 = w29*(X_1_0 + X_1_2);  
                                 const double tmp13 = w27*(X_0_0 + X_0_1);  
                                 const double tmp14 = w28*(X_1_1 + X_1_3);  
                                 const double tmp15 = w26*(X_0_2 + X_0_3);  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0 + tmp1 + tmp2 + tmp3;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp4 + tmp5 + tmp6 + tmp7;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp10 + tmp11 + tmp8 + tmp9;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp12 + tmp13 + tmp14 + tmp15;  
                             }  
                         } 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)];  
                                 EM_F[INDEX2(k,0,numEq)]+=3*X_0*w19 + 6*X_1*w20;  
                                 EM_F[INDEX2(k,1,numEq)]+=-3*X_0*w19 + 6*X_1*w20;  
                                 EM_F[INDEX2(k,2,numEq)]+=3*X_0*w19 - 6*X_1*w20;  
                                 EM_F[INDEX2(k,3,numEq)]+=-3*X_0*w19 - 6*X_1*w20;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // 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 = 6*w23*(Y_1 + Y_2);  
                                 const double tmp1 = 6*w23*(Y_0 + Y_3);  
                                 EM_F[INDEX2(k,0,numEq)]+=6*Y_0*w21 + 6*Y_3*w22 + tmp0;  
                                 EM_F[INDEX2(k,1,numEq)]+=6*Y_1*w21 + 6*Y_2*w22 + tmp1;  
                                 EM_F[INDEX2(k,2,numEq)]+=6*Y_1*w22 + 6*Y_2*w21 + tmp1;  
                                 EM_F[INDEX2(k,3,numEq)]+=6*Y_0*w22 + 6*Y_3*w21 + tmp0;  
                             }  
                         } else { // constant data  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double Y_0 = Y_p[k];  
                                 EM_F[INDEX2(k,0,numEq)]+=36*Y_0*w23;  
                                 EM_F[INDEX2(k,1,numEq)]+=36*Y_0*w23;  
                                 EM_F[INDEX2(k,2,numEq)]+=36*Y_0*w23;  
                                 EM_F[INDEX2(k,3,numEq)]+=36*Y_0*w23;  
                             }  
                         }  
                     }  
                     /* GENERATOR SNIP_PDE_SYSTEM BOTTOM */  
1943    
1944                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)  All local calculation is done on an array called `src`, with
1945                      const index_t firstNode=m_NN[0]*k1+k0;  dimensions = ext[0] * ext[1].
1946                      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  
 }  
1947    
1948  //protected  Now for MPI there is overlap to deal with. We need to share both the overlapping
1949  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  
 {  
     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;  
     const double w1 = .125*m_dx[0];  
     const double w2 = .125*m_dx[1];  
     const double w3 = .0625*m_dx[0]*m_dx[1];  
     const double w4 = .25*m_dx[0]/m_dx[1];  
     const double w5 = .25*m_dx[1]/m_dx[0];  
1950    
1951      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_NE[1]; k1+=2) {  
                 for (index_t k0=0; k0<m_NE[0]; ++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_NE[0]*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 = (A_01 + A_10)*w0;  
                                 const double tmp1 = A_00*w5;  
                                 const double tmp2 = A_01*w0;  
                                 const double tmp3 = A_10*w0;  
                                 const double tmp4 = A_11*w4;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+= tmp4 + tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+= tmp4 - tmp1 + tmp3 - tmp2;  
                                 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+= tmp2 - tmp3 - tmp4 + tmp1;  
                                 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-tmp1 - tmp4 - tmp0;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= tmp4 - tmp1 + tmp2 - tmp3;  
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+= tmp4 + tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-tmp1 + tmp0 - tmp4;  
                                 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;  
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-tmp1 + tmp0 - tmp4;  
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= tmp4 + tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= tmp4 - tmp1 + tmp2 - tmp3;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-tmp1 - tmp4 - tmp0;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= tmp2 - tmp3 - tmp4 + tmp1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= tmp4 - tmp1 + tmp3 - tmp2;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= tmp4 + tmp0 + tmp1;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process B //  
                     ///////////////  
                     if (!B.isEmpty()) {  
                         add_EM_S=true;  
                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             for (index_t m=0; m<numComp; m++) {  
                                 const double tmp0 = B_p[INDEX3(k,0,m, numEq, 2)]*w2;  
                                 const double tmp1 = B_p[INDEX3(k,1,m, numEq, 2)]*w1;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-tmp0 - tmp1;  
                                 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-tmp1 + tmp0;  
                                 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+= tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-tmp0 - tmp1;  
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-tmp1 + tmp0;  
                                 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-tmp0 - tmp1;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-tmp1 + tmp0;  
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-tmp0 - tmp1;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-tmp1 + tmp0;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= tmp0 + tmp1;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process C //  
                     ///////////////  
                     if (!C.isEmpty()) {  
                         add_EM_S=true;  
                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             for (index_t m=0; m<numComp; m++) {  
                                 const double tmp0 = C_p[INDEX3(k, m, 0, numEq, numComp)]*w2;  
                                 const double tmp1 = C_p[INDEX3(k, m, 1, numEq, numComp)]*w1;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= tmp0 - tmp1;  
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+= tmp0 - tmp1;  
                                 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+= tmp0 - tmp1;  
                                 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= tmp0 - tmp1;  
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+= tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+= tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= tmp1 - tmp0;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+= tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= tmp0 + tmp1;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= tmp0 + tmp1;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process D //  
                     ///////////////  
                     if (!D.isEmpty()) {  
                         add_EM_S=true;  
                         const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             for (index_t m=0; m<numComp; m++) {  
                                 const double tmp0 = D_p[INDEX2(k, m, numEq)]*w3;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process X //  
                     ///////////////  
                     if (!X.isEmpty()) {  
                         add_EM_F=true;  
                         const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             const double X_0 = X_p[INDEX2(k, 0, numEq)];  
                             const double X_1 = X_p[INDEX2(k, 1, numEq)];  
                             const double tmp0 = 4*X_0*w2;  
                             const double tmp1 = 4*X_1*w1;  
                             EM_F[INDEX2(k,0,numEq)]+=-tmp0 - tmp1;  
                             EM_F[INDEX2(k,1,numEq)]+= tmp0 - tmp1;  
                             EM_F[INDEX2(k,2,numEq)]+= tmp1 - tmp0;  
                             EM_F[INDEX2(k,3,numEq)]+= tmp0 + tmp1;  
                         }  
                     }  
                     ///////////////  
                     // process Y //  
                     ///////////////  
                     if (!Y.isEmpty()) {  
                         add_EM_F=true;  
                         const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             EM_F[INDEX2(k,0,numEq)]+=4*Y_p[k]*w3;  
                             EM_F[INDEX2(k,1,numEq)]+=4*Y_p[k]*w3;  
                             EM_F[INDEX2(k,2,numEq)]+=4*Y_p[k]*w3;  
                             EM_F[INDEX2(k,3,numEq)]+=4*Y_p[k]*w3;  
                         }  
                     }  
1952    
                     // add to matrix (if add_EM_S) and RHS (if add_EM_F)  
                     const index_t firstNode=m_NN[0]*k1+k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
                             firstNode, numEq, numComp);  
                 } // end k0 loop  
             } // end k1 loop  
         } // end of colouring  
     } // end of parallel region  
 }  
1953    
1954  //protected  1234567
1955  void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,  would be split into two ranks thus:
1956        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const  123(4)  (4)567     [4 being a shared element]
 {  
     const double w0 = 0.31100423396407310779*m_dx[1];  
     const double w1 = 0.022329099369260225539*m_dx[1];  
     const double w10 = 0.022329099369260225539*m_dx[0];  
     const double w11 = 0.16666666666666666667*m_dx[0];  
     const double w12 = 0.33333333333333333333*m_dx[0];  
     const double w13 = 0.39433756729740644113*m_dx[0];  
     const double w14 = 0.10566243270259355887*m_dx[0];  
     const double w15 = 0.5*m_dx[0];  
     const double w2 = 0.083333333333333333333*m_dx[1];  
     const double w3 = 0.33333333333333333333*m_dx[1];  
     const double w4 = 0.16666666666666666667*m_dx[1];  
     const double w5 = 0.39433756729740644113*m_dx[1];  
     const double w6 = 0.10566243270259355887*m_dx[1];  
     const double w7 = 0.5*m_dx[1];  
     const double w8 = 0.083333333333333333333*m_dx[0];  
     const double w9 = 0.31100423396407310779*m_dx[0];  
     const bool add_EM_S=!d.isEmpty();  
     const bool add_EM_F=!y.isEmpty();  
     rhs.requireWrite();  
 #pragma omp parallel  
     {  
         if (m_faceOffset[0] > -1) {  
             for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring  
 #pragma omp for  
                 for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     const index_t e = k1;  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  
                         if (d.actsExpanded()) {  
                             const double d_0 = d_p[0];  
                             const double d_1 = d_p[1];  
                             const double tmp0_0 = d_0 + d_1;  
                             const double tmp1_1 = d_1*w1;  
                             const double tmp4_1 = d_0*w1;  
                             const double tmp0_1 = d_0*w0;  
                             const double tmp3_1 = d_1*w0;  
                             const double tmp2_1 = tmp0_0*w2;  
                             EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp2_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp2_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;  
                         } else { /* constant data */  
                             const double d_0 = d_p[0];  
                             const double tmp1_1 = d_0*w4;  
                             const double tmp0_1 = d_0*w3;  
                             EM_S[INDEX2(0,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp1_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp0_1;  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);  
                         if (y.actsExpanded()) {  
                             const double y_0 = y_p[0];  
                             const double y_1 = y_p[1];  
                             const double tmp3_1 = w5*y_1;  
                             const double tmp2_1 = w6*y_0;  
                             const double tmp0_1 = w6*y_1;  
                             const double tmp1_1 = w5*y_0;  
                             EM_F[0]+=tmp0_1 + tmp1_1;  
                             EM_F[2]+=tmp2_1 + tmp3_1;  
                         } else { /* constant data */  
                             const double y_0 = y_p[0];  
                             const double tmp0_1 = w7*y_0;  
                             EM_F[0]+=tmp0_1;  
                             EM_F[2]+=tmp0_1;  
                         }  
                     }  
                     const index_t firstNode=m_NN[0]*k1;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 }  
             } // end colouring  
         }  
1957    
1958          if (m_faceOffset[1] > -1) {  If the radius is 2.   There will be padding elements on the outside:
             for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              
 #pragma omp for  
                 for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     const index_t e = m_faceOffset[1]+k1;  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  
                         if (d.actsExpanded()) {  
                             const double d_0 = d_p[0];  
                             const double d_1 = d_p[1];  
                             const double tmp0_0 = d_0 + d_1;  
                             const double tmp4_1 = d_1*w1;  
                             const double tmp1_1 = d_0*w1;  
                             const double tmp3_1 = d_0*w0;  
                             const double tmp0_1 = d_1*w0;  
                             const double tmp2_1 = tmp0_0*w2;  
                             EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp2_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp1_1;  
                         } else { /* constant data */  
                             const double d_0 = d_p[0];  
                             const double tmp1_1 = d_0*w4;  
                             const double tmp0_1 = d_0*w3;  
                             EM_S[INDEX2(1,1,4)]+=tmp0_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp1_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp0_1;  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);  
                         if (y.actsExpanded()) {  
                             const double y_0 = y_p[0];  
                             const double y_1 = y_p[1];  
                             const double tmp3_1 = w5*y_1;  
                             const double tmp2_1 = w6*y_0;  
                             const double tmp0_1 = w6*y_1;  
                             const double tmp1_1 = w5*y_0;  
                             EM_F[1]+=tmp0_1 + tmp1_1;  
                             EM_F[3]+=tmp2_1 + tmp3_1;  
                         } else { /* constant data */  
                             const double y_0 = y_p[0];  
                             const double tmp0_1 = w7*y_0;  
                             EM_F[1]+=tmp0_1;  
                             EM_F[3]+=tmp0_1;  
                         }  
                     }  
                     const index_t firstNode=m_NN[0]*(k1+1)-2;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 }  
             } // end colouring  
         }  
1959    
1960          if (m_faceOffset[2] > -1) {  pp123(4)  (4)567pp
             for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring  
 #pragma omp for  
                 for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     const index_t e = m_faceOffset[2]+k0;  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  
                         if (d.actsExpanded()) {  
                             const double d_0 = d_p[0];  
                             const double d_1 = d_p[1];  
                             const double tmp0_0 = d_0 + d_1;  
                             const double tmp4_1 = d_1*w9;  
                             const double tmp2_1 = d_0*w9;  
                             const double tmp0_1 = tmp0_0*w8;  
                             const double tmp1_1 = d_1*w10;  
                             const double tmp3_1 = d_0*w10;  
                             EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;  
                         } else { /* constant data */  
                             const double d_0 = d_p[0];  
                             const double tmp0_1 = d_0*w11;  
                             const double tmp1_1 = d_0*w12;  
                             EM_S[INDEX2(0,0,4)]+=tmp1_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp1_1;  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);  
                         if (y.actsExpanded()) {  
                             const double y_0 = y_p[0];  
                             const double y_1 = y_p[1];  
                             const double tmp2_1 = w13*y_1;  
                             const double tmp1_1 = w14*y_1;  
                             const double tmp3_1 = w14*y_0;  
                             const double tmp0_1 = w13*y_0;  
                             EM_F[0]+=tmp0_1 + tmp1_1;  
                             EM_F[1]+=tmp2_1 + tmp3_1;  
                         } else { /* constant data */  
                             const double y_0 = y_p[0];  
                             const double tmp0_1 = w15*y_0;  
                             EM_F[0]+=tmp0_1;  
                             EM_F[1]+=tmp0_1;  
                         }  
                     }  
                     const index_t firstNode=k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 }  
             } // end colouring  
         }  
1961    
1962          if (m_faceOffset[3] > -1) {  To ensure that 4 can be correctly computed on both ranks, values from the other rank
1963              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring  need to be known.
 #pragma omp for  
                 for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {  
                     const index_t e = m_faceOffset[3]+k0;  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  
                         if (d.actsExpanded()) {  
                             const double d_0 = d_p[0];  
                             const double d_1 = d_p[1];  
                             const double tmp0_0 = d_0 + d_1;  
                             const double tmp2_1 = d_1*w9;  
                             const double tmp4_1 = d_0*w9;  
                             const double tmp0_1 = tmp0_0*w8;  
                             const double tmp3_1 = d_1*w10;  
                             const double tmp1_1 = d_0*w10;  
                             EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;  
                         } else { /* constant data */  
                             const double d_0 = d_p[0];  
                             const double tmp0_1 = d_0*w11;  
                             const double tmp1_1 = d_0*w12;  
                             EM_S[INDEX2(2,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp1_1;  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);  
                         if (y.actsExpanded()) {  
                             const double y_0 = y_p[0];  
                             const double y_1 = y_p[1];  
                             const double tmp2_1 = w13*y_1;  
                             const double tmp1_1 = w14*y_1;  
                             const double tmp3_1 = w14*y_0;  
                             const double tmp0_1 = w13*y_0;  
                             EM_F[2]+=tmp0_1 + tmp1_1;  
                             EM_F[3]+=tmp2_1 + tmp3_1;  
                         } else { /* constant data */  
                             const double y_0 = y_p[0];  
                             const double tmp0_1 = w15*y_0;  
                             EM_F[2]+=tmp0_1;  
                             EM_F[3]+=tmp0_1;  
                         }  
                     }  
                     const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 }  
             } // end colouring  
         }  
     } // end of parallel section  
 }  
1964    
1965  //protected  pp123(4)56   23(4)567pp
 void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,  
       escript::Data& rhs, const escript::Data& d, const escript::Data& y) const  
 {  
     const double w0 = 0.25*m_dx[0];  
     const double w1 = 0.25*m_dx[1];  
     const bool add_EM_S=!d.isEmpty();  
     const bool add_EM_F=!y.isEmpty();  
     rhs.requireWrite();  
 #pragma omp parallel  
     {  
         if (m_faceOffset[0] > -1) {  
             for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring  
 #pragma omp for  
                 for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(k1);  
                         const double tmp0 = d_p[0]*w1;  
                         EM_S[INDEX2(0,0,4)]+=tmp0;  
                         EM_S[INDEX2(2,0,4)]+=tmp0;  
                         EM_S[INDEX2(0,2,4)]+=tmp0;  
                         EM_S[INDEX2(2,2,4)]+=tmp0;  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(k1);  
                         const double tmp0 = 2*w1*y_p[0];  
                         EM_F[0]+=tmp0;  
                         EM_F[2]+=tmp0;  
                     }  
                     const index_t firstNode=m_NN[0]*k1;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 }  
             } // end colouring  
         }  
1966    
1967          if (m_faceOffset[1] > -1) {  Now in our case, we wout set all the values 23456 on the left rank and send them to the
1968              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              right hand rank.
 #pragma omp for  
                 for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     const index_t e = m_faceOffset[1]+k1;  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  
                         const double tmp0 = d_p[0]*w1;  
                         EM_S[INDEX2(1,1,4)]+=tmp0;  
                         EM_S[INDEX2(3,1,4)]+=tmp0;  
                         EM_S[INDEX2(1,3,4)]+=tmp0;  
                         EM_S[INDEX2(3,3,4)]+=tmp0;  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);  
                         const double tmp0 = 2*w1*y_p[0];  
                         EM_F[1]+=tmp0;  
                         EM_F[3]+=tmp0;  
                     }  
                     const index_t firstNode=m_NN[0]*(k1+1)-2;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 }  
             } // end colouring  
         }  
1969    
1970          if (m_faceOffset[2] > -1) {  So the edges _may_ need to be shared at a distance `inset` from all boundaries.
1971              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring  inset=2*radius+1    
1972  #pragma omp for  This is to ensure that values at distance `radius` from the shared/overlapped element
1973                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {  that ripley has.
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     const index_t e = m_faceOffset[2]+k0;  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  
                         const double tmp0 = d_p[0]*w0;  
                         EM_S[INDEX2(0,0,4)]+=tmp0;  
                         EM_S[INDEX2(1,0,4)]+=tmp0;  
                         EM_S[INDEX2(0,1,4)]+=tmp0;  
                         EM_S[INDEX2(1,1,4)]+=tmp0;  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);  
                         const double tmp0 = 2*w0*y_p[0];  
                         EM_F[0]+=tmp0;  
                         EM_F[1]+=tmp0;  
                     }  
                     const index_t firstNode=k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 }  
             } // end colouring  
         }  
1974    
         if (m_faceOffset[3] > -1) {  
             for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring  
 #pragma omp for  
                 for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     const index_t e = m_faceOffset[3]+k0;  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  
                         const double tmp0 = d_p[0]*w0;  
                         EM_S[INDEX2(2,2,4)]+=tmp0;  
                         EM_S[INDEX2(3,2,4)]+=tmp0;  
                         EM_S[INDEX2(2,3,4)]+=tmp0;  
                         EM_S[INDEX2(3,3,4)]+=tmp0;  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);  
                         const double tmp0 = 2*w0*y_p[0];  
                         EM_F[2]+=tmp0;  
                         EM_F[3]+=tmp0;  
                     }  
                     const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);  
                 }  
             } // end colouring  
         }  
     } // end of parallel section  
 }  
1975    
1976  //protected  */
1977  void Rectangle::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,  escript::Data Rectangle::randomFillWorker(const escript::DataTypes::ShapeType& shape,
1978        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const         long seed, const boost::python::tuple& filter) const
1979  {  {
1980      dim_t numEq, numComp;      if (m_numDim!=2)
1981      if (!mat) {      {
1982          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());          throw RipleyException("Only 2D supported at this time.");
     } else {  
         numEq=mat->logical_row_block_size;  
         numComp=mat->logical_col_block_size;  
1983      }      }
1984      const double w0 = 0.31100423396407310779*m_dx[1];  
1985      const double w1 = 0.022329099369260225539*m_dx[1];      unsigned int radius=0;      // these are only used by gaussian
1986      const double w10 = 0.022329099369260225539*m_dx[0];      double sigma=0.5;
1987      const double w11 = 0.16666666666666666667*m_dx[0];      
1988      const double w12 = 0.33333333333333333333*m_dx[0];      unsigned int numvals=escript::DataTypes::noValues(shape);
1989      const double w13 = 0.39433756729740644113*m_dx[0];          
1990      const double w14 = 0.10566243270259355887*m_dx[0];      
1991      const double w15 = 0.5*m_dx[0];      if (len(filter)==0)
     const double w2 = 0.083333333333333333333*m_dx[1];  
     const double w3 = 0.33333333333333333333*m_dx[1];  
     const double w4 = 0.16666666666666666667*m_dx[1];  
     const double w5 = 0.39433756729740644113*m_dx[1];  
     const double w6 = 0.10566243270259355887*m_dx[1];  
     const double w7 = 0.5*m_dx[1];  
     const double w8 = 0.083333333333333333333*m_dx[0];  
     const double w9 = 0.31100423396407310779*m_dx[0];  
     const bool add_EM_S=!d.isEmpty();  
     const bool add_EM_F=!y.isEmpty();  
     rhs.requireWrite();  
 #pragma omp parallel  
1992      {      {
1993          if (m_faceOffset[0] > -1) {          // nothing special required here yet
1994              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring      }    
1995  #pragma omp for      else if (len(filter)==3)
1996                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {      {
1997                      vector<double> EM_S(4*4*numEq*numComp, 0);          boost::python::extract<string> ex(filter[0]);
1998                      vector<double> EM_F(4*numEq, 0);          if (!ex.check() || (ex()!="gaussian"))
1999                      const index_t e = k1;          {
2000                      ///////////////              throw RipleyException("Unsupported random filter");
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         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 tmp0_0 = d_0 + d_1;  
                                     const double tmp0_1 = d_0*w0;  
                                     const double tmp1_1 = d_1*w1;  
                                     const double tmp2_1 = tmp0_0*w2;  
                                     const double tmp3_1 = d_1*w0;  
                                     const double tmp4_1 = d_0*w1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_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 tmp0_1 = d_0*w3;  
                                     const double tmp1_1 = d_0*w4;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         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 tmp3_1 = w5*y_1;  
                                 const double tmp2_1 = w6*y_0;  
                                 const double tmp0_1 = w6*y_1;  
                                 const double tmp1_1 = w5*y_0;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp2_1 + tmp3_1;  
                             }  
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double y_0 = y_p[k];  
                                 const double tmp0_1 = w7*y_0;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp0_1;  
                             }  
                         }  
                     }  
                     const index_t firstNode=m_NN[0]*k1;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
                             firstNode, numEq, numComp);  
                 }  
             } // end colouring  
2001          }          }
2002            boost::python::extract<unsigned int> ex1(filter[1]);
2003          if (m_faceOffset[1] > -1) {          if (!ex1.check())
2004              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                      {
2005  #pragma omp for              throw RipleyException("Radius of gaussian filter must be a positive integer.");
                 for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {  
                     vector<double> EM_S(4*4*numEq*numComp, 0);  
                     vector<double> EM_F(4*numEq, 0);  
                     const index_t e = m_faceOffset[1]+k1;  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         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 tmp0_0 = d_0 + d_1;  
                                     const double tmp4_1 = d_1*w1;  
                                     const double tmp1_1 = d_0*w1;  
                                     const double tmp3_1 = d_0*w0;  
                                     const double tmp0_1 = d_1*w0;  
                                     const double tmp2_1 = tmp0_0*w2;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp1_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 tmp1_1 = d_0*w4;  
                                     const double tmp0_1 = d_0*w3;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         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 tmp3_1 = w5*y_1;  
                                 const double tmp2_1 = w6*y_0;  
                                 const double tmp0_1 = w6*y_1;  
                                 const double tmp1_1 = w5*y_0;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp1_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;  
                             }  
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double y_0 = y_p[k];  
                                 const double tmp0_1 = w7*y_0;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
                             }  
                         }  
                     }  
                     const index_t firstNode=m_NN[0]*(k1+1)-2;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
                             firstNode, numEq, numComp);  
                 }  
             } // end colouring  
2006          }          }
2007            radius=ex1();
2008            sigma=0.5;
2009            boost::python::extract<double> ex2(filter[2]);
2010            if (!ex2.check() || (sigma=ex2())<=0)
2011            {
2012                throw RipleyException("Sigma must be a postive floating point number.");
2013            }        
2014        }
2015        else
2016        {
2017            throw RipleyException("Unsupported random filter for Rectangle.");
2018        }
2019          
2020      
2021        
2022        size_t internal[2];
2023        internal[0]=m_NE[0]+1;      // number of points in the internal region
2024        internal[1]=m_NE[1]+1;      // that is, the ones we need smoothed versions of
2025        size_t ext[2];
2026        ext[0]=internal[0]+2*radius;        // includes points we need as input
2027        ext[1]=internal[1]+2*radius;        // for smoothing
2028        
2029        // now we check to see if the radius is acceptable
2030        // That is, would not cross multiple ranks in MPI
2031    
2032          if (m_faceOffset[2] > -1) {      if (2*radius>=internal[0]-4)
2033              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring      {
2034  #pragma omp for          throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
2035                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {      }
2036                      vector<double> EM_S(4*4*numEq*numComp, 0);      if (2*radius>=internal[1]-4)
2037                      vector<double> EM_F(4*numEq, 0);      {
2038                      const index_t e = m_faceOffset[2]+k0;          throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
2039                      ///////////////      }
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         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 tmp0_0 = d_0 + d_1;  
                                     const double tmp4_1 = d_1*w9;  
                                     const double tmp2_1 = d_0*w9;  
                                     const double tmp0_1 = tmp0_0*w8;  
                                     const double tmp1_1 = d_1*w10;  
                                     const double tmp3_1 = d_0*w10;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_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 tmp0_1 = d_0*w11;  
                                     const double tmp1_1 = d_0*w12;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         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 tmp2_1 = w13*y_1;  
                                 const double tmp1_1 = w14*y_1;  
                                 const double tmp3_1 = w14*y_0;  
                                 const double tmp0_1 = w13*y_0;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp2_1 + tmp3_1;  
                             }  
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double y_0 = y_p[k];  
                                 const double tmp0_1 = w15*y_0;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp0_1;  
                             }  
                         }  
                     }  
                     const index_t firstNode=k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
                             firstNode, numEq, numComp);  
                 }  
             } // end colouring  
         }  
2040    
2041          if (m_faceOffset[3] > -1) {      double* src=new double[ext[0]*ext[1]*numvals];
2042              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*numvals);  
2043  #pragma omp for      
                 for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {  
                     vector<double> EM_S(4*4*numEq*numComp, 0);  
                     vector<double> EM_F(4*numEq, 0);  
                     const index_t e = m_faceOffset[3]+k0;  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         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 tmp0_0 = d_0 + d_1;  
                                     const double tmp2_1 = d_1*w9;  
                                     const double tmp4_1 = d_0*w9;  
                                     const double tmp0_1 = tmp0_0*w8;  
                                     const double tmp3_1 = d_1*w10;  
                                     const double tmp1_1 = d_0*w10;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1 + tmp2_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 tmp0_1 = d_0*w11;  
                                     const double tmp1_1 = d_0*w12;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         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 tmp2_1 = w13*y_1;  
                                 const double tmp1_1 = w14*y_1;  
                                 const double tmp3_1 = w14*y_0;  
                                 const double tmp0_1 = w13*y_0;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp1_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;  
                             }  
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double y_0 = y_p[k];  
                                 const double tmp0_1 = w15*y_0;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
                             }  
                         }  
                     }  
                     const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
                             firstNode, numEq, numComp);  
                 }  
             } // end colouring  
         }  
     } // end of parallel section  
 }  
2044    
 //protected  
 void Rectangle::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,  
       escript::Data& rhs, const escript::Data& d, const escript::Data& y) const  
 {  
     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.25*m_dx[0];  
     const double w1 = 0.25*m_dx[1];  
     const bool add_EM_S=!d.isEmpty();  
     const bool add_EM_F=!y.isEmpty();  
     rhs.requireWrite();  
 #pragma omp parallel  
     {  
         if (m_faceOffset[0] > -1) {  
             for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring  
 #pragma omp for  
                 for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {  
                     vector<double> EM_S(4*4*numEq*numComp, 0);  
                     vector<double> EM_F(4*numEq, 0);  
                     ///////////////  
                     // process d //  
                     ///////////////  
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(k1);  
                         for (index_t k=0; k<numEq; k++) {  
                             for (index_t m=0; m<numComp; m++) {  
                                 const double tmp0 = d_p[INDEX2(k, m, numEq)]*w1;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(k1);  
                         for (index_t k=0; k<numEq; k++) {  
                             const double tmp0 = 2*w1*y_p[k];  
                             EM_F[INDEX2(k,0,numEq)]+=tmp0;  
                             EM_F[INDEX2(k,2,numEq)]+=tmp0;  
                         }  
                     }  
                     const index_t firstNode=m_NN[0]*k1;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
                             firstNode, numEq, numComp);  
                 }  
             } // end colouring  
         }  
2045    
2046          if (m_faceOffset[1] > -1) {  #ifdef ESYS_MPI
2047              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                  if ((internal[0]<5) || (internal[1]<5))
2048  #pragma omp for      {
2049                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {          // since the dimensions are equal for all ranks, this exception
2050                      vector<double> EM_S(4*4*numEq*numComp, 0);          // will be thrown on all ranks
2051                      vector<double> EM_F(4*numEq, 0);          throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");
2052                      const index_t e = m_faceOffset[1]+k1;      }
2053                      ///////////////      dim_t X=m_mpiInfo->rank%m_NX[0];
2054                      // process d //      dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
2055                      ///////////////  #endif      
2056                      if (add_EM_S) {  
2057                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  /*    
2058                          for (index_t k=0; k<numEq; k++) {      // if we wanted to test a repeating pattern
2059                              for (index_t m=0; m<numComp; m++) {      size_t basex=0;
2060                                  const double tmp0 = d_p[INDEX2(k, m, numEq)]*w1;      size_t basey=0;
2061                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0;  #ifdef ESYS_MPI    
2062                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0;      basex=X*m_gNE[0]/m_NX[0];
2063                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0;      basey=Y*m_gNE[1]/m_NX[1];
2064                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0;  #endif
2065                              }          
2066                          }      esysUtils::patternFillArray2D(ext[0], ext[1], src, 4, basex, basey, numvals);
2067                      }  */
2068                      ///////////////  
2069                      // process y //      
2070                      ///////////////  #ifdef ESYS_MPI  
2071                      if (add_EM_F) {      
2072                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);      BlockGrid2 grid(m_NX[0]-1, m_NX[1]-1);
2073                          for (index_t k=0; k<numEq; k++) {      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
2074                              const double tmp0 = 2*w1*y_p[k];                  // there is an overlap of two points both of which need to have "radius" points on either side.
2075                              EM_F[INDEX2(k,1,numEq)]+=tmp0;      
2076                              EM_F[INDEX2(k,3,numEq)]+=tmp0;      size_t xmidlen=ext[0]-2*inset;      // how wide is the x-dimension between the two insets
2077                          }      size_t ymidlen=ext[1]-2*inset;      
2078                      }      
2079                      const index_t firstNode=m_NN[0]*(k1+1)-2;      Block2 block(ext[0], ext[1], inset, xmidlen, ymidlen, numvals);
2080                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
2081                              firstNode, numEq, numComp);      MPI_Request reqs[40];               // a non-tight upper bound on how many we need
2082                  }      MPI_Status stats[40];
2083              } // end colouring      short rused=0;
2084          }      
2085        messvec incoms;
2086        messvec outcoms;  
2087        
2088        grid.generateInNeighbours(X, Y, incoms);
2089        grid.generateOutNeighbours(X, Y, outcoms);
2090        
2091        block.copyAllToBuffer(src);  
2092        
2093        int comserr=0;    
2094        for (size_t i=0;i<incoms.size();++i)
2095        {
2096            message& m=incoms[i];
2097            comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
2098            block.setUsed(m.destbuffid);
2099        
2100        }
2101    
2102          if (m_faceOffset[2] > -1) {      for (size_t i=0;i<outcoms.size();++i)
2103              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring      {
2104  #pragma omp for          message& m=outcoms[i];
2105                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {          comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
2106                      vector<double> EM_S(4*4*numEq*numComp, 0);      }    
2107                      vector<double> EM_F(4*numEq, 0);      
2108                      const index_t e = m_faceOffset[2]+k0;      if (!comserr)
2109                      ///////////////      {    
2110                      // process d //          comserr=MPI_Waitall(rused, reqs, stats);    
2111                      ///////////////      }
                     if (add_EM_S) {  
                         const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             for (index_t m=0; m<numComp; m++) {  
                                 const double tmp0 = d_p[INDEX2(k, m, numEq)]*w0;  
                                 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0;  
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0;  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process y //  
                     ///////////////  
                     if (add_EM_F) {  
                         const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);  
                         for (index_t k=0; k<numEq; k++) {  
                             const double tmp0 = 2*w0*y_p[k];  
                             EM_F[INDEX2(k,0,numEq)]+=tmp0;  
                             EM_F[INDEX2(k,1,numEq)]+=tmp0;  
                         }  
                     }  
                     const index_t firstNode=k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
                             firstNode, numEq, numComp);  
                 }  
             } // end colouring  
         }  
2112    
2113          if (m_faceOffset[3] > -1) {      if (comserr)
2114              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring      {
2115  #pragma omp for          // Yes this is throwing an exception as a result of an MPI error.
2116                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {          // and no we don't inform the other ranks that we are doing this.
2117                      vector<double> EM_S(4*4*numEq*numComp, 0);          // however, we have no reason to believe coms work at this point anyway
2118                      vector<double> EM_F(4*numEq, 0);          throw RipleyException("Error in coms for randomFill");      
2119                      const index_t e = m_faceOffset[3]+k0;      }
2120                      ///////////////      
2121                      // process d //      block.copyUsedFromBuffer(src);    
2122                      ///////////////      
2123                      if (add_EM_S) {  #endif    
2124                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);      
2125                          for (index_t k=0; k<numEq; k++) {      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
2126                              for (index_t m=0; m<numComp; m++) {      {
2127                                  const double tmp0 = d_p[INDEX2(k, m, numEq)]*w0;        
2128                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0;          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
2129                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0;          escript::Data resdat(0, shape, fs , true);
2130                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0;          // don't need to check for exwrite because we just made it
2131                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0;          escript::DataVector& dv=resdat.getExpandedVectorReference();
2132                              }          
2133                          }          
2134                      }          // now we need to copy values over
2135                      ///////////////          for (size_t y=0;y<(internal[1]);++y)    
2136                      // process y //          {
2137                      ///////////////              for (size_t x=0;x<(internal[0]);++x)
2138                      if (add_EM_F) {              {
2139                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                  for (unsigned int i=0;i<numvals;++i)
2140                          for (index_t k=0; k<numEq; k++) {                  {
2141                              const double tmp0 = 2*w0*y_p[k];                      dv[i+(x+y*(internal[0]))*numvals]=src[i+(x+y*ext[0])*numvals];
                             EM_F[INDEX2(k,2,numEq)]+=tmp0;  
                             EM_F[INDEX2(k,3,numEq)]+=tmp0;  
                         }  
                     }  
                     const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;  
                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,  
                             firstNode, numEq, numComp);  
2142                  }                  }
2143              } // end colouring              }
2144          }          }
2145      } // end of parallel section          delete[] src;
2146            return resdat;      
2147        }
2148        else                // filter enabled      
2149        {    
2150            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
2151            escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
2152            // don't need to check for exwrite because we just made it
2153            escript::DataVector& dv=resdat.getExpandedVectorReference();
2154            double* convolution=get2DGauss(radius, sigma);
2155            for (size_t y=0;y<(internal[1]);++y)    
2156            {
2157                for (size_t x=0;x<(internal[0]);++x)
2158                {    
2159                    dv[x+y*(internal[0])]=Convolve2D(convolution, src, x+radius, y+radius, radius, ext[0]);
2160                    
2161                }
2162            }
2163            delete[] convolution;
2164            delete[] src;
2165            return resdat;
2166        }
2167    }
2168    
2169    int Rectangle::findNode(const double *coords) const {
2170        const int NOT_MINE = -1;
2171        //is the found element even owned by this rank
2172        // (inside owned or shared elements but will map to an owned element)
2173        for (int dim = 0; dim < m_numDim; dim++) {
2174            double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
2175                    - m_dx[dim]/2.; //allows for point outside mapping onto node
2176            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
2177                    + m_dx[dim]/2.;
2178            if (min > coords[dim] || max < coords[dim]) {
2179                return NOT_MINE;
2180            }
2181        }
2182        // get distance from origin
2183        double x = coords[0] - m_origin[0];
2184        double y = coords[1] - m_origin[1];
2185        // distance in elements
2186        int ex = (int) floor(x / m_dx[0]);
2187        int ey = (int) floor(y / m_dx[1]);
2188        // set the min distance high enough to be outside the element plus a bit
2189        int closest = NOT_MINE;
2190        double minDist = 1;
2191        for (int dim = 0; dim < m_numDim; dim++) {
2192            minDist += m_dx[dim]*m_dx[dim];
2193        }
2194        //find the closest node
2195        for (int dx = 0; dx < 1; dx++) {
2196            double xdist = (x - (ex + dx)*m_dx[0]);
2197            for (int dy = 0; dy < 1; dy++) {
2198                double ydist = (y - (ey + dy)*m_dx[1]);
2199                double total = xdist*xdist + ydist*ydist;
2200                if (total < minDist) {
2201                    closest = INDEX2(ex+dx-m_offset[0], ey+dy-m_offset[1], m_NE[0] + 1);
2202                    minDist = total;
2203                }
2204            }
2205        }
2206        //if this happens, we've let a dirac point slip through, which is awful
2207        if (closest == NOT_MINE) {
2208            throw RipleyException("Unable to map appropriate dirac point to a node,"
2209                    " implementation problem in Rectangle::findNode()");
2210        }
2211        return closest;
2212    }
2213    
2214    void Rectangle::setAssembler(std::string type, std::map<std::string,
2215            escript::Data> constants) {
2216        if (type.compare("WaveAssembler") == 0) {
2217            delete assembler;
2218            assembler = new WaveAssembler2D(this, m_dx, m_NX, m_NE, m_NN, constants);
2219        } else { //else ifs would go before this for other types
2220            throw RipleyException("Ripley::Rectangle does not support the"
2221                                    " requested assembler");
2222        }
2223  }  }
2224    
   
2225  } // end of namespace ripley  } // end of namespace ripley
2226    

Legend:
Removed from v.4370  
changed lines
  Added in v.4705

  ViewVC Help
Powered by ViewVC 1.1.26