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

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

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

branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp revision 3762 by caltinay, Tue Jan 10 00:01:45 2012 UTC trunk/ripley/src/Rectangle.cpp revision 4705 by jfenwick, Fri Feb 21 02:36:15 2014 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2011 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development 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  extern "C" {  #include <paso/SystemMatrix.h>
19  #include "paso/SystemMatrix.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
27    #include <netcdfcpp.h>
28    #endif
29    
30  #if USE_SILO  #if USE_SILO
31  #include <silo.h>  #include <silo.h>
# Line 26  extern "C" { Line 37  extern "C" {
37  #include <iomanip>  #include <iomanip>
38    
39  using namespace std;  using namespace std;
40    using esysUtils::FileWriter;
41    
42  namespace ripley {  namespace ripley {
43    
44  Rectangle::Rectangle(int n0, int n1, double l0, double l1, int d0, int d1) :  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
45      RipleyDomain(2),                       double y1, int d0, int d1,
46      m_gNE0(n0),                       const std::vector<double>& points,
47      m_gNE1(n1),                       const std::vector<int>& tags,
48      m_l0(l0),                       const simap_t& tagnamestonums) :
49      m_l1(l1),      RipleyDomain(2)
50      m_NX(d0),  {
51      m_NY(d1)      // ignore subdivision parameters for serial run
52  {      if (m_mpiInfo->size == 1) {
53            d0=1;
54            d1=1;
55        }
56    
57        bool warn=false;
58        // if number of subdivisions is non-positive, try to subdivide by the same
59        // ratio as the number of elements
60        if (d0<=0 && d1<=0) {
61            warn=true;
62            d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
63            d1=m_mpiInfo->size/d0;
64            if (d0*d1 != m_mpiInfo->size) {
65                // ratios not the same so subdivide side with more elements only
66                if (n0>n1) {
67                    d0=0;
68                    d1=1;
69                } else {
70                    d0=1;
71                    d1=0;
72                }
73            }
74        }
75        if (d0<=0) {
76            // d1 is preset, determine d0 - throw further down if result is no good
77            d0=m_mpiInfo->size/d1;
78        } else if (d1<=0) {
79            // d0 is preset, determine d1 - throw further down if result is no good
80            d1=m_mpiInfo->size/d0;
81        }
82    
83      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
84      // among number of ranks      // among number of ranks
85      if (m_NX*m_NY != m_mpiInfo->size)      if (d0*d1 != m_mpiInfo->size)
86          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
87    
88      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)      if (warn) {
89          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
90                << d1 << "). This may not be optimal!" << endl;
91        }
92    
93        double l0 = x1-x0;
94        double l1 = y1-y0;
95        m_dx[0] = l0/n0;
96        m_dx[1] = l1/n1;
97    
98        if ((n0+1)%d0 > 0) {
99            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
100            l0=m_dx[0]*n0;
101            cout << "Warning: Adjusted number of elements and length. N0="
102                << n0 << ", l0=" << l0 << endl;
103        }
104        if ((n1+1)%d1 > 0) {
105            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
106            l1=m_dx[1]*n1;
107            cout << "Warning: Adjusted number of elements and length. N1="
108                << n1 << ", l1=" << l1 << endl;
109        }
110    
111      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))      if ((d0 > 1 && (n0+1)/d0<2) || (d1 > 1 && (n1+1)/d1<2))
112          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
113    
114      // local number of elements (including overlap)      m_gNE[0] = n0;
115      m_NE0 = (m_NX>1 ? (n0+1)/m_NX : n0);      m_gNE[1] = n1;
116      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)      m_origin[0] = x0;
117          m_NE0++;      m_origin[1] = y0;
118      m_NE1 = (m_NY>1 ? (n1+1)/m_NY : n1);      m_length[0] = l0;
119      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)      m_length[1] = l1;
120          m_NE1++;      m_NX[0] = d0;
121        m_NX[1] = d1;
122    
123        // local number of elements (with and without overlap)
124        m_NE[0] = m_ownNE[0] = (d0>1 ? (n0+1)/d0 : n0);
125        if (m_mpiInfo->rank%d0>0 && m_mpiInfo->rank%d0<d0-1)
126            m_NE[0]++;
127        else if (d0>1 && m_mpiInfo->rank%d0==d0-1)
128            m_ownNE[0]--;
129    
130        m_NE[1] = m_ownNE[1] = (d1>1 ? (n1+1)/d1 : n1);
131        if (m_mpiInfo->rank/d0>0 && m_mpiInfo->rank/d0<d1-1)
132            m_NE[1]++;
133        else if (d1>1 && m_mpiInfo->rank/d0==d1-1)
134            m_ownNE[1]--;
135    
136      // local number of nodes      // local number of nodes
137      m_N0 = m_NE0+1;      m_NN[0] = m_NE[0]+1;
138      m_N1 = m_NE1+1;      m_NN[1] = m_NE[1]+1;
139    
140      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
141      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);      m_offset[0] = (n0+1)/d0*(m_mpiInfo->rank%d0);
142      if (m_offset0 > 0)      if (m_offset[0] > 0)
143          m_offset0--;          m_offset[0]--;
144      m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);      m_offset[1] = (n1+1)/d1*(m_mpiInfo->rank/d0);
145      if (m_offset1 > 0)      if (m_offset[1] > 0)
146          m_offset1--;          m_offset[1]--;
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);
161        Paso_Connector_free(m_connector);
162        delete assembler;
163  }  }
164    
165  string Rectangle::getDescription() const  string Rectangle::getDescription() const
# Line 87  bool Rectangle::operator==(const Abstrac Line 172  bool Rectangle::operator==(const Abstrac
172      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
173      if (o) {      if (o) {
174          return (RipleyDomain::operator==(other) &&          return (RipleyDomain::operator==(other) &&
175                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1                  m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1]
176                  && m_l0==o->m_l0 && m_l1==o->m_l1                  && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1]
177                  && m_NX==o->m_NX && m_NY==o->m_NY);                  && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1]
178                    && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1]);
179      }      }
180    
181      return false;      return false;
182  }  }
183    
184    void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
185                const ReaderParameters& params) const
186    {
187    #ifdef USE_NETCDF
188        // check destination function space
189        int myN0, myN1;
190        if (out.getFunctionSpace().getTypeCode() == Nodes) {
191            myN0 = m_NN[0];
192            myN1 = m_NN[1];
193        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
194                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
195            myN0 = m_NE[0];
196            myN1 = m_NE[1];
197        } else
198            throw RipleyException("readNcGrid(): invalid function space for output data object");
199    
200        if (params.first.size() != 2)
201            throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
202    
203        if (params.numValues.size() != 2)
204            throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
205    
206        if (params.multiplier.size() != 2)
207            throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");
208        for (size_t i=0; i<params.multiplier.size(); i++)
209            if (params.multiplier[i]<1)
210                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
215        NcFile f(filename.c_str(), NcFile::ReadOnly);
216        if (!f.is_valid())
217            throw RipleyException("readNcGrid(): cannot open file");
218    
219        NcVar* var = f.get_var(varname.c_str());
220        if (!var)
221            throw RipleyException("readNcGrid(): invalid variable");
222    
223        // TODO: rank>0 data support
224        const int numComp = out.getDataPointSize();
225        if (numComp > 1)
226            throw RipleyException("readNcGrid(): only scalar data supported");
227    
228        const int dims = var->num_dims();
229        boost::scoped_array<long> edges(var->edges());
230    
231        // is this a slice of the data object (dims!=2)?
232        // note the expected ordering of edges (as in numpy: y,x)
233        if ( (dims==2 && (params.numValues[1] > edges[0] || params.numValues[0] > edges[1]))
234                || (dims==1 && params.numValues[1]>1) ) {
235            throw RipleyException("readNcGrid(): not enough data in file");
236        }
237    
238        // check if this rank contributes anything
239        if (params.first[0] >= m_offset[0]+myN0 ||
240                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;
244    
245        // now determine how much this rank has to write
246    
247        // first coordinates in data object to write to
248        const int first0 = max(0, params.first[0]-m_offset[0]);
249        const int first1 = max(0, params.first[1]-m_offset[1]);
250        // indices to first value in file (not accounting for reverse yet)
251        int idx0 = max(0, m_offset[0]-params.first[0]);
252        int idx1 = max(0, m_offset[1]-params.first[1]);
253        // number of values to read
254        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
255        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);
264        if (dims==2) {
265            var->set_cur(idx1, idx0);
266            var->get(&values[0], num1, num0);
267        } else {
268            var->set_cur(idx0);
269            var->get(&values[0], num0);
270        }
271    
272        const int dpp = out.getNumDataPointsPerSample();
273        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++) {
282    #pragma omp parallel for
283            for (index_t x=0; x<num0; x++) {
284                const int baseIndex = first0+x*params.multiplier[0]
285                                      +(first1+y*params.multiplier[1])*myN0;
286                const int srcIndex = (y0+y_mult*y)*num0+(x0+x_mult*x);
287                if (!isnan(values[srcIndex])) {
288                    for (index_t m1=0; m1<params.multiplier[1]; m1++) {
289                        for (index_t m0=0; m0<params.multiplier[0]; m0++) {
290                            const int dataIndex = baseIndex+m0+m1*myN0;
291                            double* dest = out.getSampleDataRW(dataIndex);
292                            for (index_t q=0; q<dpp; q++) {
293                                *dest++ = values[srcIndex];
294                            }
295                        }
296                    }
297                }
298            }
299        }
300    #else
301        throw RipleyException("readNcGrid(): not compiled with netCDF support");
302    #endif
303    }
304    
305    void Rectangle::readBinaryGrid(escript::Data& out, string filename,
306                                   const ReaderParameters& params) const
307    {
308        // 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
330        int myN0, myN1;
331        if (out.getFunctionSpace().getTypeCode() == Nodes) {
332            myN0 = m_NN[0];
333            myN1 = m_NN[1];
334        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
335                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
336            myN0 = m_NE[0];
337            myN1 = m_NE[1];
338        } else
339            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
340    
341        // check file existence and size
342        ifstream f(filename.c_str(), ifstream::binary);
343        if (f.fail()) {
344            throw RipleyException("readBinaryGrid(): cannot open file");
345        }
346        f.seekg(0, ios::end);
347        const int numComp = out.getDataPointSize();
348        const int filesize = f.tellg();
349        const int reqsize = params.numValues[0]*params.numValues[1]*numComp*sizeof(ValueType);
350        if (filesize < reqsize) {
351            f.close();
352            throw RipleyException("readBinaryGrid(): not enough data in file");
353        }
354    
355        // check if this rank contributes anything
356        if (params.first[0] >= m_offset[0]+myN0 ||
357                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();
361            return;
362        }
363    
364        // now determine how much this rank has to write
365    
366        // first coordinates in data object to write to
367        const int first0 = max(0, params.first[0]-m_offset[0]);
368        const int first1 = max(0, params.first[1]-m_offset[1]);
369        // indices to first value in file
370        const int idx0 = max(0, m_offset[0]-params.first[0]);
371        const int idx1 = max(0, m_offset[1]-params.first[1]);
372        // number of values to read
373        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
374        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
375    
376        out.requireWrite();
377        vector<ValueType> values(num0*numComp);
378        const int dpp = out.getNumDataPointsPerSample();
379    
380        for (int y=0; y<num1; y++) {
381            const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]);
382            f.seekg(fileofs*sizeof(ValueType));
383            f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
384            for (int x=0; x<num0; x++) {
385                const int baseIndex = first0+x*params.multiplier[0]
386                                        +(first1+y*params.multiplier[1])*myN0;
387                for (int m1=0; m1<params.multiplier[1]; m1++) {
388                    for (int m0=0; m0<params.multiplier[0]; m0++) {
389                        const int dataIndex = baseIndex+m0+m1*myN0;
390                        double* dest = out.getSampleDataRW(dataIndex);
391                        for (int c=0; c<numComp; c++) {
392                            ValueType val = values[x*numComp+c];
393    
394                            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                        }
405                    }
406                }
407            }
408        }
409    
410        f.close();
411    }
412    
413    void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,
414                                    int byteOrder, int dataType) const
415    {
416        // the mapping is not universally correct but should work on our
417        // supported platforms
418        switch (dataType) {
419            case DATATYPE_INT32:
420                writeBinaryGridImpl<int>(in, filename, byteOrder);
421                break;
422            case DATATYPE_FLOAT32:
423                writeBinaryGridImpl<float>(in, filename, byteOrder);
424                break;
425            case DATATYPE_FLOAT64:
426                writeBinaryGridImpl<double>(in, filename, byteOrder);
427                break;
428            default:
429                throw RipleyException("writeBinaryGrid(): invalid or unsupported datatype");
430        }
431    }
432    
433    template<typename ValueType>
434    void Rectangle::writeBinaryGridImpl(const escript::Data& in,
435                                        const string& filename, int byteOrder) const
436    {
437        // check function space and determine number of points
438        int myN0, myN1;
439        int totalN0, totalN1;
440        if (in.getFunctionSpace().getTypeCode() == Nodes) {
441            myN0 = m_NN[0];
442            myN1 = m_NN[1];
443            totalN0 = m_gNE[0]+1;
444            totalN1 = m_gNE[1]+1;
445        } else if (in.getFunctionSpace().getTypeCode() == Elements ||
446                    in.getFunctionSpace().getTypeCode() == ReducedElements) {
447            myN0 = m_NE[0];
448            myN1 = m_NE[1];
449            totalN0 = m_gNE[0];
450            totalN1 = m_gNE[1];
451        } else
452            throw RipleyException("writeBinaryGrid(): invalid function space of data object");
453    
454        const int numComp = in.getDataPointSize();
455        const int dpp = in.getNumDataPointsPerSample();
456    
457        if (numComp > 1 || dpp > 1)
458            throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
459    
460        const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;
461    
462        // from here on we know that each sample consists of one value
463        FileWriter fw;
464        fw.openFile(filename, fileSize);
465        MPIBarrier();
466    
467        for (index_t y=0; y<myN1; y++) {
468            const int fileofs = (m_offset[0]+(m_offset[1]+y)*totalN0)*sizeof(ValueType);
469            ostringstream oss;
470    
471            for (index_t x=0; x<myN0; x++) {
472                const double* sample = in.getSampleDataRO(y*myN0+x);
473                ValueType fvalue = static_cast<ValueType>(*sample);
474                if (byteOrder == BYTEORDER_NATIVE) {
475                    oss.write((char*)&fvalue, sizeof(fvalue));
476                } else {
477                    char* value = reinterpret_cast<char*>(&fvalue);
478                    oss.write(byte_swap32(value), sizeof(fvalue));
479                }
480            }
481            fw.writeAt(oss, fileofs);
482        }
483        fw.close();
484    }
485    
486  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
487  {  {
488  #if USE_SILO  #if USE_SILO
# Line 103  void Rectangle::dump(const string& fileN Line 491  void Rectangle::dump(const string& fileN
491          fn+=".silo";          fn+=".silo";
492      }      }
493    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
494      int driver=DB_HDF5;          int driver=DB_HDF5;    
     string siloPath;  
495      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
496        const char* blockDirFmt = "/block%04d";
497    
498  #ifdef ESYS_MPI  #ifdef ESYS_MPI
499      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
500        const int NUM_SILO_FILES = 1;
501  #endif  #endif
502    
503      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 126  void Rectangle::dump(const string& fileN Line 513  void Rectangle::dump(const string& fileN
513                          PMPIO_DefaultClose, (void*)&driver);                          PMPIO_DefaultClose, (void*)&driver);
514          }          }
515          if (baton) {          if (baton) {
516              char str[64];              char siloPath[64];
517              snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));              snprintf(siloPath, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
518              siloPath = str;              dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath);
             dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());  
519          }          }
520  #endif  #endif
521      } else {      } else {
# Line 141  void Rectangle::dump(const string& fileN Line 527  void Rectangle::dump(const string& fileN
527              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
528                      getDescription().c_str(), driver);                      getDescription().c_str(), driver);
529          }          }
530            char siloPath[64];
531            snprintf(siloPath, 64, blockDirFmt, 0);
532            DBMkDir(dbfile, siloPath);
533            DBSetDir(dbfile, siloPath);
534      }      }
535    
536      if (!dbfile)      if (!dbfile)
# Line 155  void Rectangle::dump(const string& fileN Line 545  void Rectangle::dump(const string& fileN
545      }      }
546      */      */
547    
548      boost::scoped_ptr<double> x(new double[m_N0]);      boost::scoped_ptr<double> x(new double[m_NN[0]]);
549      boost::scoped_ptr<double> y(new double[m_N1]);      boost::scoped_ptr<double> y(new double[m_NN[1]]);
550      double* coords[2] = { x.get(), y.get() };      double* coords[2] = { x.get(), y.get() };
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
551  #pragma omp parallel  #pragma omp parallel
552      {      {
553  #pragma omp for nowait  #pragma omp for nowait
554          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
555              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=getLocalCoordinate(i0, 0);
556          }          }
557  #pragma omp for nowait  #pragma omp for nowait
558          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
559              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=getLocalCoordinate(i1, 1);
560          }          }
561      }      }
562      IndexVector dims = getNumNodesPerDim();      int* dims = const_cast<int*>(getNumNodesPerDim());
563    
564      // write mesh      // write mesh
565      DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,
566              DB_COLLINEAR, NULL);              DB_COLLINEAR, NULL);
567    
568      // write node ids      // write node ids
569      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 2,      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 2,
570              NULL, 0, DB_INT, DB_NODECENT, NULL);              NULL, 0, DB_INT, DB_NODECENT, NULL);
571    
572      // write element ids      // write element ids
573      dims = getNumElementsPerDim();      dims = const_cast<int*>(getNumElementsPerDim());
574      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
575              &dims[0], 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);              dims, 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);
576    
577      // rank 0 writes multimesh and multivar      // rank 0 writes multimesh and multivar
578      if (m_mpiInfo->rank == 0) {      if (m_mpiInfo->rank == 0) {
# Line 233  void Rectangle::dump(const string& fileN Line 621  void Rectangle::dump(const string& fileN
621      }      }
622    
623  #else // USE_SILO  #else // USE_SILO
624      throw RipleyException("dump(): no Silo support");      throw RipleyException("dump: no Silo support");
625  #endif  #endif
626  }  }
627    
# Line 241  const int* Rectangle::borrowSampleRefere Line 629  const int* Rectangle::borrowSampleRefere
629  {  {
630      switch (fsType) {      switch (fsType) {
631          case Nodes:          case Nodes:
632          case ReducedNodes: //FIXME: reduced          case ReducedNodes: // FIXME: reduced
633              return &m_nodeId[0];              return &m_nodeId[0];
634          case DegreesOfFreedom:          case DegreesOfFreedom:
635          case ReducedDegreesOfFreedom: //FIXME: reduced          case ReducedDegreesOfFreedom: // FIXME: reduced
636              return &m_dofId[0];              return &m_dofId[0];
637          case Elements:          case Elements:
638          case ReducedElements:          case ReducedElements:
# Line 252  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      }      }
648    
649      stringstream msg;      stringstream msg;
650      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
651      throw RipleyException(msg.str());      throw RipleyException(msg.str());
652  }  }
653    
# Line 269  bool Rectangle::ownSample(int fsType, in Line 658  bool Rectangle::ownSample(int fsType, in
658    
659      switch (fsType) {      switch (fsType) {
660          case Nodes:          case Nodes:
661          case ReducedNodes: //FIXME: reduced          case ReducedNodes: // FIXME: reduced
662              return (m_dofMap[id] < getNumDOF());              return (m_dofMap[id] < getNumDOF());
663          case DegreesOfFreedom:          case DegreesOfFreedom:
664          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
# Line 277  bool Rectangle::ownSample(int fsType, in Line 666  bool Rectangle::ownSample(int fsType, in
666          case Elements:          case Elements:
667          case ReducedElements:          case ReducedElements:
668              // check ownership of element's bottom left node              // check ownership of element's bottom left node
669              return (m_dofMap[id%m_NE0+m_N0*(id/m_NE0)] < getNumDOF());              return (m_dofMap[id%m_NE[0]+m_NN[0]*(id/m_NE[0])] < getNumDOF());
670          case FaceElements:          case FaceElements:
671          case ReducedFaceElements:          case ReducedFaceElements:
672              {              {
673                  // check ownership of face element's first node                  // determine which face the sample belongs to before
674                  const IndexVector faces = getNumFacesPerBoundary();                  // checking ownership of corresponding element's first node
675                  dim_t n=0;                  dim_t n=0;
676                  for (size_t i=0; i<faces.size(); i++) {                  for (size_t i=0; i<4; i++) {
677                      n+=faces[i];                      n+=m_faceCount[i];
678                      if (id<n) {                      if (id<n) {
679                          index_t k;                          index_t k;
680                          if (i==1)                          if (i==1)
681                              k=m_N0-1;                              k=m_NN[0]-2;
682                          else if (i==3)                          else if (i==3)
683                              k=m_N0*(m_N1-1);                              k=m_NN[0]*(m_NN[1]-2);
684                          else                          else
685                              k=0;                              k=0;
686                          // determine whether to move right or up                          // determine whether to move right or up
687                          const index_t delta=(i/2==0 ? m_N0 : 1);                          const index_t delta=(i/2==0 ? m_NN[0] : 1);
688                          return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());                          return (m_dofMap[k+(id-n+m_faceCount[i])*delta] < getNumDOF());
689                      }                      }
690                  }                  }
691                  return false;                  return false;
# Line 306  bool Rectangle::ownSample(int fsType, in Line 695  bool Rectangle::ownSample(int fsType, in
695      }      }
696    
697      stringstream msg;      stringstream msg;
698      msg << "ownSample() not implemented for "      msg << "ownSample: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
699      throw RipleyException(msg.str());      throw RipleyException(msg.str());
700  }  }
701    
 void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const  
 {  
     escript::Data& in = *const_cast<escript::Data*>(&cIn);  
     const dim_t numComp = in.getDataPointSize();  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     const double cx0 = -1./h0;  
     const double cx1 = -.78867513459481288225/h0;  
     const double cx2 = -.5/h0;  
     const double cx3 = -.21132486540518711775/h0;  
     const double cx4 = .21132486540518711775/h0;  
     const double cx5 = .5/h0;  
     const double cx6 = .78867513459481288225/h0;  
     const double cx7 = 1./h0;  
     const double cy0 = -1./h1;  
     const double cy1 = -.78867513459481288225/h1;  
     const double cy2 = -.5/h1;  
     const double cy3 = -.21132486540518711775/h1;  
     const double cy4 = .21132486540518711775/h1;  
     const double cy5 = .5/h1;  
     const double cy6 = .78867513459481288225/h1;  
     const double cy7 = 1./h1;  
   
     if (out.getFunctionSpace().getTypeCode() == Elements) {  
         out.requireWrite();  
         /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */  
 #pragma omp parallel for  
         for (index_t k1=0; k1 < m_NE1; ++k1) {  
             for (index_t k0=0; k0 < m_NE0; ++k0) {  
                 const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));  
                 const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));  
                 const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  
                 const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  
                 double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));  
                 for (index_t i=0; i < numComp; ++i) {  
                     o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;  
                     o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;  
                     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,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;  
                     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,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;  
                     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,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;  
                 } /* end of component loop i */  
             } /* end of k0 loop */  
         } /* end of k1 loop */  
         /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */  
     } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {  
         out.requireWrite();  
         /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */  
 #pragma omp parallel for  
         for (index_t k1=0; k1 < m_NE1; ++k1) {  
             for (index_t k0=0; k0 < m_NE0; ++k0) {  
                 const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));  
                 const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));  
                 const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  
                 const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  
                 double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));  
                 for (index_t i=0; i < numComp; ++i) {  
                     o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);  
                     o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);  
                 } /* end of component loop i */  
             } /* end of k0 loop */  
         } /* end of k1 loop */  
         /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */  
     } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {  
         out.requireWrite();  
 #pragma omp parallel  
         {  
             /*** GENERATOR SNIP_GRAD_FACES TOP */  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1=0; k1 < m_NE1; ++k1) {  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;  
                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;  
                         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,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;  
                     } /* end of component loop i */  
                 } /* end of k1 loop */  
             } /* end of face 0 */  
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1=0; k1 < m_NE1; ++k1) {  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;  
                         o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;  
                         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,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;  
                     } /* end of component loop i */  
                 } /* end of k1 loop */  
             } /* end of face 1 */  
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;  
                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;  
                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;  
                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of face 2 */  
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;  
                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;  
                         o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;  
                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of face 3 */  
             /* GENERATOR SNIP_GRAD_FACES BOTTOM */  
         } // end of parallel section  
     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
         out.requireWrite();  
 #pragma omp parallel  
         {  
             /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1=0; k1 < m_NE1; ++k1) {  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);  
                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;  
                     } /* end of component loop i */  
                 } /* end of k1 loop */  
             } /* end of face 0 */  
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1=0; k1 < m_NE1; ++k1) {  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);  
                         o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;  
                     } /* end of component loop i */  
                 } /* end of k1 loop */  
             } /* end of face 1 */  
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;  
                         o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of face 2 */  
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;  
                         o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of face 3 */  
             /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */  
         } // end of parallel section  
     } else {  
         stringstream msg;  
         msg << "setToGradient() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
 }  
   
 void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  
 {  
     escript::Data& in = *const_cast<escript::Data*>(&arg);  
     const dim_t numComp = in.getDataPointSize();  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     if (arg.getFunctionSpace().getTypeCode() == Elements) {  
         const double w_0 = h0*h1/4.;  
 #pragma omp parallel  
         {  
             vector<double> int_local(numComp, 0);  
 #pragma omp for nowait  
             for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         const register double f_2 = f[INDEX2(i,2,numComp)];  
                         const register double f_3 = f[INDEX2(i,3,numComp)];  
                         int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k1 loop */  
   
 #pragma omp critical  
             for (index_t i=0; i<numComp; i++)  
                 integrals[i]+=int_local[i];  
         } // end of parallel section  
     } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
         const double w_0 = h0*h1;  
 #pragma omp parallel  
         {  
             vector<double> int_local(numComp, 0);  
 #pragma omp for nowait  
             for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*w_0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k1 loop */  
   
 #pragma omp critical  
             for (index_t i=0; i<numComp; i++)  
                 integrals[i]+=int_local[i];  
         } // end of parallel section  
     } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
         const double w_0 = h0/2.;  
         const double w_1 = h1/2.;  
 #pragma omp parallel  
         {  
             vector<double> int_local(numComp, 0);  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         int_local[i]+=(f_0+f_1)*w_1;  
                     }  /* end of component loop i */  
                 } /* end of k1 loop */  
             }  
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         int_local[i]+=(f_0+f_1)*w_1;  
                     }  /* end of component loop i */  
                 } /* end of k1 loop */  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         int_local[i]+=(f_0+f_1)*w_0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         int_local[i]+=(f_0+f_1)*w_0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             }  
   
 #pragma omp critical  
             for (index_t i=0; i<numComp; i++)  
                 integrals[i]+=int_local[i];  
         } // end of parallel section  
     } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
 #pragma omp parallel  
         {  
             vector<double> int_local(numComp, 0);  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*h1;  
                     }  /* end of component loop i */  
                 } /* end of k1 loop */  
             }  
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*h1;  
                     }  /* end of component loop i */  
                 } /* end of k1 loop */  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*h0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*h0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             }  
   
 #pragma omp critical  
             for (index_t i=0; i<numComp; i++)  
                 integrals[i]+=int_local[i];  
         } // end of parallel section  
     } else {  
         stringstream msg;  
         msg << "setToIntegrals() not implemented for "  
             << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
 }  
   
702  void Rectangle::setToNormal(escript::Data& out) const  void Rectangle::setToNormal(escript::Data& out) const
703  {  {
704      if (out.getFunctionSpace().getTypeCode() == FaceElements) {      if (out.getFunctionSpace().getTypeCode() == FaceElements) {
# Line 688  void Rectangle::setToNormal(escript::Dat Line 707  void Rectangle::setToNormal(escript::Dat
707          {          {
708              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
709  #pragma omp for nowait  #pragma omp for nowait
710                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
711                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
712                      // set vector at two quadrature points                      // set vector at two quadrature points
713                      *o++ = -1.;                      *o++ = -1.;
# Line 700  void Rectangle::setToNormal(escript::Dat Line 719  void Rectangle::setToNormal(escript::Dat
719    
720              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
721  #pragma omp for nowait  #pragma omp for nowait
722                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
723                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
724                      // set vector at two quadrature points                      // set vector at two quadrature points
725                      *o++ = 1.;                      *o++ = 1.;
# Line 712  void Rectangle::setToNormal(escript::Dat Line 731  void Rectangle::setToNormal(escript::Dat
731    
732              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
733  #pragma omp for nowait  #pragma omp for nowait
734                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
735                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
736                      // set vector at two quadrature points                      // set vector at two quadrature points
737                      *o++ = 0.;                      *o++ = 0.;
# Line 724  void Rectangle::setToNormal(escript::Dat Line 743  void Rectangle::setToNormal(escript::Dat
743    
744              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
745  #pragma omp for nowait  #pragma omp for nowait
746                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
747                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
748                      // set vector at two quadrature points                      // set vector at two quadrature points
749                      *o++ = 0.;                      *o++ = 0.;
# Line 740  void Rectangle::setToNormal(escript::Dat Line 759  void Rectangle::setToNormal(escript::Dat
759          {          {
760              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
761  #pragma omp for nowait  #pragma omp for nowait
762                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
763                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
764                      *o++ = -1.;                      *o++ = -1.;
765                      *o = 0.;                      *o = 0.;
# Line 749  void Rectangle::setToNormal(escript::Dat Line 768  void Rectangle::setToNormal(escript::Dat
768    
769              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
770  #pragma omp for nowait  #pragma omp for nowait
771                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
772                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
773                      *o++ = 1.;                      *o++ = 1.;
774                      *o = 0.;                      *o = 0.;
# Line 758  void Rectangle::setToNormal(escript::Dat Line 777  void Rectangle::setToNormal(escript::Dat
777    
778              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
779  #pragma omp for nowait  #pragma omp for nowait
780                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
781                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
782                      *o++ = 0.;                      *o++ = 0.;
783                      *o = -1.;                      *o = -1.;
# Line 767  void Rectangle::setToNormal(escript::Dat Line 786  void Rectangle::setToNormal(escript::Dat
786    
787              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
788  #pragma omp for nowait  #pragma omp for nowait
789                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
790                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
791                      *o++ = 0.;                      *o++ = 0.;
792                      *o = 1.;                      *o = 1.;
# Line 777  void Rectangle::setToNormal(escript::Dat Line 796  void Rectangle::setToNormal(escript::Dat
796    
797      } else {      } else {
798          stringstream msg;          stringstream msg;
799          msg << "setToNormal() not implemented for "          msg << "setToNormal: invalid function space type "
800              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());              << out.getFunctionSpace().getTypeCode();
801          throw RipleyException(msg.str());          throw RipleyException(msg.str());
802      }      }
803  }  }
# Line 789  void Rectangle::setToSize(escript::Data& Line 808  void Rectangle::setToSize(escript::Data&
808              || out.getFunctionSpace().getTypeCode() == ReducedElements) {              || out.getFunctionSpace().getTypeCode() == ReducedElements) {
809          out.requireWrite();          out.requireWrite();
810          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
811          const double hSize=getFirstCoordAndSpacing(0).second;          const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]);
         const double vSize=getFirstCoordAndSpacing(1).second;  
         const double size=min(hSize,vSize);  
812  #pragma omp parallel for  #pragma omp parallel for
813          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
814              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 801  void Rectangle::setToSize(escript::Data& Line 818  void Rectangle::setToSize(escript::Data&
818              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
819          out.requireWrite();          out.requireWrite();
820          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
         const double hSize=getFirstCoordAndSpacing(0).second;  
         const double vSize=getFirstCoordAndSpacing(1).second;  
821  #pragma omp parallel  #pragma omp parallel
822          {          {
823              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
824  #pragma omp for nowait  #pragma omp for nowait
825                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
826                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
827                      fill(o, o+numQuad, vSize);                      fill(o, o+numQuad, m_dx[1]);
828                  }                  }
829              }              }
830    
831              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
832  #pragma omp for nowait  #pragma omp for nowait
833                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
834                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
835                      fill(o, o+numQuad, vSize);                      fill(o, o+numQuad, m_dx[1]);
836                  }                  }
837              }              }
838    
839              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
840  #pragma omp for nowait  #pragma omp for nowait
841                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
842                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
843                      fill(o, o+numQuad, hSize);                      fill(o, o+numQuad, m_dx[0]);
844                  }                  }
845              }              }
846    
847              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
848  #pragma omp for nowait  #pragma omp for nowait
849                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
850                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
851                      fill(o, o+numQuad, hSize);                      fill(o, o+numQuad, m_dx[0]);
852                  }                  }
853              }              }
854          } // end of parallel section          } // end of parallel section
855    
856      } else {      } else {
857          stringstream msg;          stringstream msg;
858          msg << "setToSize() not implemented for "          msg << "setToSize: invalid function space type "
859              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());              << out.getFunctionSpace().getTypeCode();
860          throw RipleyException(msg.str());          throw RipleyException(msg.str());
861      }      }
862  }  }
863    
 Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  
                                                 bool reducedColOrder) const  
 {  
     if (reducedRowOrder || reducedColOrder)  
         throw RipleyException("getPattern() not implemented for reduced order");  
   
     return m_pattern;  
 }  
   
864  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
865  {  {
866      RipleyDomain::Print_Mesh_Info(full);      RipleyDomain::Print_Mesh_Info(full);
# Line 862  void Rectangle::Print_Mesh_Info(const bo Line 868  void Rectangle::Print_Mesh_Info(const bo
868          cout << "     Id  Coordinates" << endl;          cout << "     Id  Coordinates" << endl;
869          cout.precision(15);          cout.precision(15);
870          cout.setf(ios::scientific, ios::floatfield);          cout.setf(ios::scientific, ios::floatfield);
         pair<double,double> xdx = getFirstCoordAndSpacing(0);  
         pair<double,double> ydy = getFirstCoordAndSpacing(1);  
871          for (index_t i=0; i < getNumNodes(); i++) {          for (index_t i=0; i < getNumNodes(); i++) {
872              cout << "  " << setw(5) << m_nodeId[i]              cout << "  " << setw(5) << m_nodeId[i]
873                  << "  " << xdx.first+(i%m_N0)*xdx.second                  << "  " << getLocalCoordinate(i%m_NN[0], 0)
874                  << "  " << ydy.first+(i/m_N0)*ydy.second << endl;                  << "  " << getLocalCoordinate(i/m_NN[0], 1) << endl;
875          }          }
876      }      }
877  }  }
878    
 IndexVector Rectangle::getNumNodesPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_N0);  
     ret.push_back(m_N1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumElementsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NE0);  
     ret.push_back(m_NE1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumFacesPerBoundary() const  
 {  
     IndexVector ret(4, 0);  
     //left  
     if (m_offset0==0)  
         ret[0]=m_NE1;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         ret[1]=m_NE1;  
     //bottom  
     if (m_offset1==0)  
         ret[2]=m_NE0;  
     //top  
     if (m_mpiInfo->rank/m_NX==m_NY-1)  
         ret[3]=m_NE0;  
     return ret;  
 }  
   
 pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  
 {  
     if (dim==0) {  
         return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);  
     } else if (dim==1) {  
         return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);  
     }  
     throw RipleyException("getFirstCoordAndSpacing(): invalid argument");  
 }  
   
 //protected  
 dim_t Rectangle::getNumDOF() const  
 {  
     return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;  
 }  
   
 //protected  
 dim_t Rectangle::getNumFaceElements() const  
 {  
     const IndexVector faces = getNumFacesPerBoundary();  
     dim_t n=0;  
     for (size_t i=0; i<faces.size(); i++)  
         n+=faces[i];  
     return n;  
 }  
879    
880  //protected  //protected
881  void Rectangle::assembleCoordinates(escript::Data& arg) const  void Rectangle::assembleCoordinates(escript::Data& arg) const
# Line 942  void Rectangle::assembleCoordinates(escr Line 887  void Rectangle::assembleCoordinates(escr
887      if (!numSamplesEqual(&x, 1, getNumNodes()))      if (!numSamplesEqual(&x, 1, getNumNodes()))
888          throw RipleyException("setToX: Illegal number of samples in Data object");          throw RipleyException("setToX: Illegal number of samples in Data object");
889    
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
890      arg.requireWrite();      arg.requireWrite();
891  #pragma omp parallel for  #pragma omp parallel for
892      for (dim_t i1 = 0; i1 < m_N1; i1++) {      for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
893          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
894              double* point = arg.getSampleDataRW(i0+m_N0*i1);              double* point = arg.getSampleDataRW(i0+m_NN[0]*i1);
895              point[0] = xdx.first+i0*xdx.second;              point[0] = getLocalCoordinate(i0, 0);
896              point[1] = ydy.first+i1*ydy.second;              point[1] = getLocalCoordinate(i1, 1);
897          }          }
898      }      }
899  }  }
900    
901  //protected  //protected
902    void Rectangle::assembleGradient(escript::Data& out, const escript::Data& in) const
903    {
904        const dim_t numComp = in.getDataPointSize();
905        const double cx0 = .21132486540518711775/m_dx[0];
906        const double cx1 = .78867513459481288225/m_dx[0];
907        const double cx2 = 1./m_dx[0];
908        const double cy0 = .21132486540518711775/m_dx[1];
909        const double cy1 = .78867513459481288225/m_dx[1];
910        const double cy2 = 1./m_dx[1];
911    
912        if (out.getFunctionSpace().getTypeCode() == Elements) {
913            out.requireWrite();
914    #pragma omp parallel
915            {
916                vector<double> f_00(numComp);
917                vector<double> f_01(numComp);
918                vector<double> f_10(numComp);
919                vector<double> f_11(numComp);
920    #pragma omp for
921                for (index_t k1=0; k1 < m_NE[1]; ++k1) {
922                    for (index_t k0=0; k0 < m_NE[0]; ++k0) {
923                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
924                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
925                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
926                        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]));
928                        for (index_t i=0; i < numComp; ++i) {
929                            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_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
931                            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_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
933                            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_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
935                            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_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
937                        } // end of component loop i
938                    } // end of k0 loop
939                } // end of k1 loop
940            } // end of parallel section
941        } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
942            out.requireWrite();
943    #pragma omp parallel
944            {
945                vector<double> f_00(numComp);
946                vector<double> f_01(numComp);
947                vector<double> f_10(numComp);
948                vector<double> f_11(numComp);
949    #pragma omp for
950                for (index_t k1=0; k1 < m_NE[1]; ++k1) {
951                    for (index_t k0=0; k0 < m_NE[0]; ++k0) {
952                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
953                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
954                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
955                        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]));
957                        for (index_t i=0; i < numComp; ++i) {
958                            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)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
960                        } // end of component loop i
961                    } // end of k0 loop
962                } // end of k1 loop
963            } // end of parallel section
964        } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
965            out.requireWrite();
966    #pragma omp parallel
967            {
968                vector<double> f_00(numComp);
969                vector<double> f_01(numComp);
970                vector<double> f_10(numComp);
971                vector<double> f_11(numComp);
972                if (m_faceOffset[0] > -1) {
973    #pragma omp for nowait
974                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
975                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
976                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
977                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
978                        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);
980                        for (index_t i=0; i < numComp; ++i) {
981                            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_01[i]-f_00[i])*cy2;
983                            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_01[i]-f_00[i])*cy2;
985                        } // end of component loop i
986                    } // end of k1 loop
987                } // end of face 0
988                if (m_faceOffset[1] > -1) {
989    #pragma omp for nowait
990                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
991                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
992                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
993                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
994                        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);
996                        for (index_t i=0; i < numComp; ++i) {
997                            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_11[i]-f_10[i])*cy2;
999                            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_11[i]-f_10[i])*cy2;
1001                        } // end of component loop i
1002                    } // end of k1 loop
1003                } // end of face 1
1004                if (m_faceOffset[2] > -1) {
1005    #pragma omp for nowait
1006                    for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1007                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1008                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
1009                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1010                        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);
1012                        for (index_t i=0; i < numComp; ++i) {
1013                            o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx2;
1014                            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_10[i]-f_00[i])*cx2;
1016                            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
1018                    } // end of k0 loop
1019                } // end of face 2
1020                if (m_faceOffset[3] > -1) {
1021    #pragma omp for nowait
1022                    for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1023                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1024                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1025                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1026                        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);
1028                        for (index_t i=0; i < numComp; ++i) {
1029                            o[INDEX3(i,0,0,numComp,2)] = (f_11[i]-f_01[i])*cx2;
1030                            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_11[i]-f_01[i])*cx2;
1032                            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
1034                    } // end of k0 loop
1035                } // end of face 3
1036            } // end of parallel section
1037    
1038        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1039            out.requireWrite();
1040    #pragma omp parallel
1041            {
1042                vector<double> f_00(numComp);
1043                vector<double> f_01(numComp);
1044                vector<double> f_10(numComp);
1045                vector<double> f_11(numComp);
1046                if (m_faceOffset[0] > -1) {
1047    #pragma omp for nowait
1048                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1049                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1050                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1051                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
1052                        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);
1054                        for (index_t i=0; i < numComp; ++i) {
1055                            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_01[i]-f_00[i])*cy2;
1057                        } // end of component loop i
1058                    } // end of k1 loop
1059                } // end of face 0
1060                if (m_faceOffset[1] > -1) {
1061    #pragma omp for nowait
1062                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1063                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
1064                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
1065                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1066                        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);
1068                        for (index_t i=0; i < numComp; ++i) {
1069                            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_11[i]-f_10[i])*cy2;
1071                        } // end of component loop i
1072                    } // end of k1 loop
1073                } // end of face 1
1074                if (m_faceOffset[2] > -1) {
1075    #pragma omp for nowait
1076                    for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1077                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1078                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
1079                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1080                        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);
1082                        for (index_t i=0; i < numComp; ++i) {
1083                            o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx2;
1084                            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
1086                    } // end of k0 loop
1087                } // end of face 2
1088                if (m_faceOffset[3] > -1) {
1089    #pragma omp for nowait
1090                    for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1091                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1092                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1093                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1094                        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);
1096                        for (index_t i=0; i < numComp; ++i) {
1097                            o[INDEX3(i,0,0,numComp,2)] = (f_11[i]-f_01[i])*cx2;
1098                            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
1100                    } // end of k0 loop
1101                } // end of face 3
1102            } // end of parallel section
1103        }
1104    }
1105    
1106    //protected
1107    void Rectangle::assembleIntegrate(vector<double>& integrals,
1108                                      const escript::Data& arg) const
1109    {
1110        const dim_t numComp = arg.getDataPointSize();
1111        const index_t left = (m_offset[0]==0 ? 0 : 1);
1112        const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1113        const int fs=arg.getFunctionSpace().getTypeCode();
1114        if (fs == Elements && arg.actsExpanded()) {
1115    #pragma omp parallel
1116            {
1117                vector<double> int_local(numComp, 0);
1118                const double w = m_dx[0]*m_dx[1]/4.;
1119    #pragma omp for nowait
1120                for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1121                    for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1122                        const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1123                        for (index_t i=0; i < numComp; ++i) {
1124                            const double f0 = f[INDEX2(i,0,numComp)];
1125                            const double f1 = f[INDEX2(i,1,numComp)];
1126                            const double f2 = f[INDEX2(i,2,numComp)];
1127                            const double f3 = f[INDEX2(i,3,numComp)];
1128                            int_local[i]+=(f0+f1+f2+f3)*w;
1129                        }  // end of component loop i
1130                    } // end of k0 loop
1131                } // end of k1 loop
1132    #pragma omp critical
1133                for (index_t i=0; i<numComp; i++)
1134                    integrals[i]+=int_local[i];
1135            } // end of parallel section
1136    
1137        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1138            const double w = m_dx[0]*m_dx[1];
1139    #pragma omp parallel
1140            {
1141                vector<double> int_local(numComp, 0);
1142    #pragma omp for nowait
1143                for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1144                    for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1145                        const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1146                        for (index_t i=0; i < numComp; ++i) {
1147                            int_local[i]+=f[i]*w;
1148                        }
1149                    }
1150                }
1151    #pragma omp critical
1152                for (index_t i=0; i<numComp; i++)
1153                    integrals[i]+=int_local[i];
1154            } // end of parallel section
1155    
1156        } else if (fs == FaceElements && arg.actsExpanded()) {
1157    #pragma omp parallel
1158            {
1159                vector<double> int_local(numComp, 0);
1160                const double w0 = m_dx[0]/2.;
1161                const double w1 = m_dx[1]/2.;
1162                if (m_faceOffset[0] > -1) {
1163    #pragma omp for nowait
1164                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1165                        const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1166                        for (index_t i=0; i < numComp; ++i) {
1167                            const double f0 = f[INDEX2(i,0,numComp)];
1168                            const double f1 = f[INDEX2(i,1,numComp)];
1169                            int_local[i]+=(f0+f1)*w1;
1170                        }  // end of component loop i
1171                    } // end of k1 loop
1172                }
1173    
1174                if (m_faceOffset[1] > -1) {
1175    #pragma omp for nowait
1176                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1177                        const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1178                        for (index_t i=0; i < numComp; ++i) {
1179                            const double f0 = f[INDEX2(i,0,numComp)];
1180                            const double f1 = f[INDEX2(i,1,numComp)];
1181                            int_local[i]+=(f0+f1)*w1;
1182                        }  // end of component loop i
1183                    } // end of k1 loop
1184                }
1185    
1186                if (m_faceOffset[2] > -1) {
1187    #pragma omp for nowait
1188                    for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1189                        const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1190                        for (index_t i=0; i < numComp; ++i) {
1191                            const double f0 = f[INDEX2(i,0,numComp)];
1192                            const double f1 = f[INDEX2(i,1,numComp)];
1193                            int_local[i]+=(f0+f1)*w0;
1194                        }  // end of component loop i
1195                    } // end of k0 loop
1196                }
1197    
1198                if (m_faceOffset[3] > -1) {
1199    #pragma omp for nowait
1200                    for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1201                        const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1202                        for (index_t i=0; i < numComp; ++i) {
1203                            const double f0 = f[INDEX2(i,0,numComp)];
1204                            const double f1 = f[INDEX2(i,1,numComp)];
1205                            int_local[i]+=(f0+f1)*w0;
1206                        }  // end of component loop i
1207                    } // end of k0 loop
1208                }
1209    #pragma omp critical
1210                for (index_t i=0; i<numComp; i++)
1211                    integrals[i]+=int_local[i];
1212            } // end of parallel section
1213    
1214        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1215    #pragma omp parallel
1216            {
1217                vector<double> int_local(numComp, 0);
1218                if (m_faceOffset[0] > -1) {
1219    #pragma omp for nowait
1220                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1221                        const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1222                        for (index_t i=0; i < numComp; ++i) {
1223                            int_local[i]+=f[i]*m_dx[1];
1224                        }
1225                    }
1226                }
1227    
1228                if (m_faceOffset[1] > -1) {
1229    #pragma omp for nowait
1230                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1231                        const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1232                        for (index_t i=0; i < numComp; ++i) {
1233                            int_local[i]+=f[i]*m_dx[1];
1234                        }
1235                    }
1236                }
1237    
1238                if (m_faceOffset[2] > -1) {
1239    #pragma omp for nowait
1240                    for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1241                        const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1242                        for (index_t i=0; i < numComp; ++i) {
1243                            int_local[i]+=f[i]*m_dx[0];
1244                        }
1245                    }
1246                }
1247    
1248                if (m_faceOffset[3] > -1) {
1249    #pragma omp for nowait
1250                    for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1251                        const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1252                        for (index_t i=0; i < numComp; ++i) {
1253                            int_local[i]+=f[i]*m_dx[0];
1254                        }
1255                    }
1256                }
1257    
1258    #pragma omp critical
1259                for (index_t i=0; i<numComp; i++)
1260                    integrals[i]+=int_local[i];
1261            } // end of parallel section
1262        } // function space selector
1263    }
1264    
1265    //protected
1266  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
1267  {  {
1268      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1269      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1270      const int x=node%nDOF0;      const int x=node%nDOF0;
1271      const int y=node/nDOF0;      const int y=node/nDOF0;
1272      dim_t num=0;      dim_t num=0;
# Line 984  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();
1298    
1299      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1300      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1301      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1302      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1303  #pragma omp parallel for  #pragma omp parallel for
1304      for (index_t i=0; i<nDOF1; i++) {      for (index_t i=0; i<nDOF1; i++) {
1305          for (index_t j=0; j<nDOF0; j++) {          for (index_t j=0; j<nDOF0; j++) {
1306              const index_t n=j+left+(i+bottom)*m_N0;              const index_t n=j+left+(i+bottom)*m_NN[0];
1307              const double* src=in.getSampleDataRO(n);              const double* src=in.getSampleDataRO(n);
1308              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
1309          }          }
# Line 1004  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 1022  void Rectangle::dofToNodes(escript::Data Line 1330  void Rectangle::dofToNodes(escript::Data
1330                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1331          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1332      }      }
1333        Paso_Coupler_free(coupler);
1334  }  }
1335    
1336  //private  //private
1337  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
1338  {  {
1339      // identifiers are ordered from left to right, bottom to top globablly.      // degrees of freedom are numbered from left to right, bottom to top in
1340        // each rank, continuing on the next rank (ranks also go left-right,
1341        // bottom-top).
1342        // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
1343        // helps when writing out data rank after rank.
1344    
1345      // build node distribution vector first.      // build node distribution vector first.
1346      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
1347        // constant for all ranks in this implementation
1348      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1349      const dim_t numDOF=getNumDOF();      const dim_t numDOF=getNumDOF();
1350      for (dim_t k=1; k<m_mpiInfo->size; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
# Line 1040  void Rectangle::populateSampleIds() Line 1354  void Rectangle::populateSampleIds()
1354      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
1355      m_dofId.resize(numDOF);      m_dofId.resize(numDOF);
1356      m_elementId.resize(getNumElements());      m_elementId.resize(getNumElements());
1357    
1358        // populate face element counts
1359        //left
1360        if (m_offset[0]==0)
1361            m_faceCount[0]=m_NE[1];
1362        else
1363            m_faceCount[0]=0;
1364        //right
1365        if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
1366            m_faceCount[1]=m_NE[1];
1367        else
1368            m_faceCount[1]=0;
1369        //bottom
1370        if (m_offset[1]==0)
1371            m_faceCount[2]=m_NE[0];
1372        else
1373            m_faceCount[2]=0;
1374        //top
1375        if (m_mpiInfo->rank/m_NX[0]==m_NX[1]-1)
1376            m_faceCount[3]=m_NE[0];
1377        else
1378            m_faceCount[3]=0;
1379    
1380      m_faceId.resize(getNumFaceElements());      m_faceId.resize(getNumFaceElements());
1381    
1382        const index_t left = (m_offset[0]==0 ? 0 : 1);
1383        const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1384        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1385        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1386    
1387    #define globalNodeId(x,y) \
1388        ((m_offset[0]+x)/nDOF0)*nDOF0*nDOF1+(m_offset[0]+x)%nDOF0 \
1389        + ((m_offset[1]+y)/nDOF1)*nDOF0*nDOF1*m_NX[0]+((m_offset[1]+y)%nDOF1)*nDOF0
1390    
1391        // set corner id's outside the parallel region
1392        m_nodeId[0] = globalNodeId(0, 0);
1393        m_nodeId[m_NN[0]-1] = globalNodeId(m_NN[0]-1, 0);
1394        m_nodeId[m_NN[0]*(m_NN[1]-1)] = globalNodeId(0, m_NN[1]-1);
1395        m_nodeId[m_NN[0]*m_NN[1]-1] = globalNodeId(m_NN[0]-1,m_NN[1]-1);
1396    #undef globalNodeId
1397    
1398  #pragma omp parallel  #pragma omp parallel
1399      {      {
1400          // nodes          // populate degrees of freedom and own nodes (identical id)
1401  #pragma omp for nowait  #pragma omp for nowait
1402          for (dim_t i1=0; i1<m_N1; i1++) {          for (dim_t i=0; i<nDOF1; i++) {
1403              for (dim_t i0=0; i0<m_N0; i0++) {              for (dim_t j=0; j<nDOF0; j++) {
1404                  m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;                  const index_t nodeIdx=j+left+(i+bottom)*m_NN[0];
1405                    const index_t dofIdx=j+i*nDOF0;
1406                    m_dofId[dofIdx] = m_nodeId[nodeIdx]
1407                        = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
1408              }              }
1409          }          }
1410    
1411          // degrees of freedom          // populate the rest of the nodes (shared with other ranks)
1412            if (m_faceCount[0]==0) { // left column
1413  #pragma omp for nowait  #pragma omp for nowait
1414          for (dim_t k=0; k<numDOF; k++)              for (dim_t i=0; i<nDOF1; i++) {
1415              m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;                  const index_t nodeIdx=(i+bottom)*m_NN[0];
1416                    const index_t dofId=(i+1)*nDOF0-1;
1417                    m_nodeId[nodeIdx]
1418                        = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
1419                }
1420            }
1421            if (m_faceCount[1]==0) { // right column
1422    #pragma omp for nowait
1423                for (dim_t i=0; i<nDOF1; i++) {
1424                    const index_t nodeIdx=(i+bottom+1)*m_NN[0]-1;
1425                    const index_t dofId=i*nDOF0;
1426                    m_nodeId[nodeIdx]
1427                        = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
1428                }
1429            }
1430            if (m_faceCount[2]==0) { // bottom row
1431    #pragma omp for nowait
1432                for (dim_t i=0; i<nDOF0; i++) {
1433                    const index_t nodeIdx=i+left;
1434                    const index_t dofId=nDOF0*(nDOF1-1)+i;
1435                    m_nodeId[nodeIdx]
1436                        = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
1437                }
1438            }
1439            if (m_faceCount[3]==0) { // top row
1440    #pragma omp for nowait
1441                for (dim_t i=0; i<nDOF0; i++) {
1442                    const index_t nodeIdx=m_NN[0]*(m_NN[1]-1)+i+left;
1443                    const index_t dofId=i;
1444                    m_nodeId[nodeIdx]
1445                        = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
1446                }
1447            }
1448    
1449          // elements          // populate element id's
1450  #pragma omp for nowait  #pragma omp for nowait
1451          for (dim_t i1=0; i1<m_NE1; i1++) {          for (dim_t i1=0; i1<m_NE[1]; i1++) {
1452              for (dim_t i0=0; i0<m_NE0; i0++) {              for (dim_t i0=0; i0<m_NE[0]; i0++) {
1453                  m_elementId[i0+i1*m_NE0]=(m_offset1+i1)*m_gNE0+m_offset0+i0;                  m_elementId[i0+i1*m_NE[0]]=(m_offset[1]+i1)*m_gNE[0]+m_offset[0]+i0;
1454              }              }
1455          }          }
1456    
# Line 1078  void Rectangle::populateSampleIds() Line 1467  void Rectangle::populateSampleIds()
1467      updateTagsInUse(Elements);      updateTagsInUse(Elements);
1468    
1469      // generate face offset vector and set face tags      // generate face offset vector and set face tags
     const IndexVector facesPerEdge = getNumFacesPerBoundary();  
1470      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
1471      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
1472      m_faceOffset.assign(facesPerEdge.size(), -1);      m_faceOffset.assign(4, -1);
1473      m_faceTags.clear();      m_faceTags.clear();
1474      index_t offset=0;      index_t offset=0;
1475      for (size_t i=0; i<facesPerEdge.size(); i++) {      for (size_t i=0; i<4; i++) {
1476          if (facesPerEdge[i]>0) {          if (m_faceCount[i]>0) {
1477              m_faceOffset[i]=offset;              m_faceOffset[i]=offset;
1478              offset+=facesPerEdge[i];              offset+=m_faceCount[i];
1479              m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);              m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
1480          }          }
1481      }      }
1482      setTagMap("left", LEFT);      setTagMap("left", LEFT);
# Line 1101  void Rectangle::populateSampleIds() Line 1489  void Rectangle::populateSampleIds()
1489  //private  //private
1490  void Rectangle::createPattern()  void Rectangle::createPattern()
1491  {  {
1492      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1493      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1494      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1495      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1496    
1497      // populate node->DOF mapping with own degrees of freedom.      // populate node->DOF mapping with own degrees of freedom.
1498      // The rest is assigned in the loop further down      // The rest is assigned in the loop further down
1499      m_dofMap.assign(getNumNodes(), 0);      m_dofMap.assign(getNumNodes(), 0);
1500  #pragma omp parallel for  #pragma omp parallel for
1501      for (index_t i=bottom; i<m_N1; i++) {      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1502          for (index_t j=left; j<m_N0; j++) {          for (index_t j=left; j<left+nDOF0; j++) {
1503              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;              m_dofMap[i*m_NN[0]+j]=(i-bottom)*nDOF0+j-left;
1504          }          }
1505      }      }
1506    
# Line 1125  void Rectangle::createPattern() Line 1513  void Rectangle::createPattern()
1513      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
1514      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
1515      int numShared=0;      int numShared=0;
1516      const int x=m_mpiInfo->rank%m_NX;      const int x=m_mpiInfo->rank%m_NX[0];
1517      const int y=m_mpiInfo->rank/m_NX;      const int y=m_mpiInfo->rank/m_NX[0];
1518      for (int i1=-1; i1<2; i1++) {      for (int i1=-1; i1<2; i1++) {
1519          for (int i0=-1; i0<2; i0++) {          for (int i0=-1; i0<2; i0++) {
1520              // skip this rank              // skip this rank
# Line 1135  void Rectangle::createPattern() Line 1523  void Rectangle::createPattern()
1523              // location of neighbour rank              // location of neighbour rank
1524              const int nx=x+i0;              const int nx=x+i0;
1525              const int ny=y+i1;              const int ny=y+i1;
1526              if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {              if (nx>=0 && ny>=0 && nx<m_NX[0] && ny<m_NX[1]) {
1527                  neighbour.push_back(ny*m_NX+nx);                  neighbour.push_back(ny*m_NX[0]+nx);
1528                  if (i0==0) {                  if (i0==0) {
1529                      // sharing top or bottom edge                      // sharing top or bottom edge
1530                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1531                      const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);                      const int firstNode=(i1==-1 ? left : m_NN[0]*(m_NN[1]-1)+left);
1532                      offsetInShared.push_back(offsetInShared.back()+nDOF0);                      offsetInShared.push_back(offsetInShared.back()+nDOF0);
1533                      for (dim_t i=0; i<nDOF0; i++, numShared++) {                      for (dim_t i=0; i<nDOF0; i++, numShared++) {
1534                          sendShared.push_back(firstDOF+i);                          sendShared.push_back(firstDOF+i);
# Line 1155  void Rectangle::createPattern() Line 1543  void Rectangle::createPattern()
1543                  } else if (i1==0) {                  } else if (i1==0) {
1544                      // sharing left or right edge                      // sharing left or right edge
1545                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1546                      const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);                      const int firstNode=(i0==-1 ? bottom*m_NN[0] : (bottom+1)*m_NN[0]-1);
1547                      offsetInShared.push_back(offsetInShared.back()+nDOF1);                      offsetInShared.push_back(offsetInShared.back()+nDOF1);
1548                      for (dim_t i=0; i<nDOF1; i++, numShared++) {                      for (dim_t i=0; i<nDOF1; i++, numShared++) {
1549                          sendShared.push_back(firstDOF+i*nDOF0);                          sendShared.push_back(firstDOF+i*nDOF0);
# Line 1165  void Rectangle::createPattern() Line 1553  void Rectangle::createPattern()
1553                          colIndices[firstDOF+i*nDOF0].push_back(numShared);                          colIndices[firstDOF+i*nDOF0].push_back(numShared);
1554                          if (i<nDOF1-1)                          if (i<nDOF1-1)
1555                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1556                          m_dofMap[firstNode+i*m_N0]=numDOF+numShared;                          m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
1557                      }                      }
1558                  } else {                  } else {
1559                      // sharing a node                      // sharing a node
1560                      const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);                      const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1561                      const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);                      const int node=(i0+1)/2*(m_NN[0]-1)+(i1+1)/2*m_NN[0]*(m_NN[1]-1);
1562                      offsetInShared.push_back(offsetInShared.back()+1);                      offsetInShared.push_back(offsetInShared.back()+1);
1563                      sendShared.push_back(dof);                      sendShared.push_back(dof);
1564                      recvShared.push_back(numDOF+numShared);                      recvShared.push_back(numDOF+numShared);
# Line 1272  void Rectangle::createPattern() Line 1660  void Rectangle::createPattern()
1660      Paso_Pattern_free(rowPattern);      Paso_Pattern_free(rowPattern);
1661  }  }
1662    
1663    //private
1664    void Rectangle::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
1665             const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
1666             bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
1667    {
1668        IndexVector rowIndex;
1669        rowIndex.push_back(m_dofMap[firstNode]);
1670        rowIndex.push_back(m_dofMap[firstNode+1]);
1671        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]]);
1672        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]+1]);
1673        if (addF) {
1674            double *F_p=F.getSampleDataRW(0);
1675            for (index_t i=0; i<rowIndex.size(); i++) {
1676                if (rowIndex[i]<getNumDOF()) {
1677                    for (index_t eq=0; eq<nEq; eq++) {
1678                        F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
1679                    }
1680                }
1681            }
1682        }
1683        if (addS) {
1684            addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
1685        }
1686    }
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) {
1695          out.requireWrite();          out.requireWrite();
1696          /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */          const double c0 = 0.25;
1697          const double c0 = .25;  #pragma omp parallel
1698  #pragma omp parallel for          {
1699          for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_00(numComp);
1700              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_01(numComp);
1701                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1702                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_11(numComp);
1703                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1704                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1705                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1706                  for (index_t i=0; i < numComp; ++i) {                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
1707                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
1708                  } /* end of component loop i */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
1709              } /* end of k0 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
1710          } /* end of k1 loop */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1711          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */                      for (index_t i=0; i < numComp; ++i) {
1712                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1713                        } /* end of component loop i */
1714                    } /* end of k0 loop */
1715                } /* end of k1 loop */
1716            } /* end of parallel section */
1717      } else {      } else {
1718          out.requireWrite();          out.requireWrite();
1719          /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */          const double c0 = 0.16666666666666666667;
1720          const double c0 = .16666666666666666667;          const double c1 = 0.044658198738520451079;
1721          const double c1 = .044658198738520451079;          const double c2 = 0.62200846792814621559;
1722          const double c2 = .62200846792814621559;  #pragma omp parallel
1723  #pragma omp parallel for          {
1724          for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_00(numComp);
1725              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_01(numComp);
1726                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1727                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_11(numComp);
1728                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1729                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1730                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1731                  for (index_t i=0; i < numComp; ++i) {                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
1732                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
1733                      o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
1734                      o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
1735                      o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1736                  } /* end of component loop i */                      for (index_t i=0; i < numComp; ++i) {
1737              } /* end of k0 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1738          } /* end of k1 loop */                          o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1739          /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */                          o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1740                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1741                        } /* end of component loop i */
1742                    } /* end of k0 loop */
1743                } /* end of k1 loop */
1744            } /* end of parallel section */
1745      }      }
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 = .5;  
1756  #pragma omp parallel  #pragma omp parallel
1757          {          {
1758              /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */              vector<double> f_00(numComp);
1759                vector<double> f_01(numComp);
1760                vector<double> f_10(numComp);
1761                vector<double> f_11(numComp);
1762              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1763  #pragma omp for nowait  #pragma omp for nowait
1764                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1765                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1766                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
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 */
1773              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1774  #pragma omp for nowait  #pragma omp for nowait
1775                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1776                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1777                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
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 */
1784              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1785  #pragma omp for nowait  #pragma omp for nowait
1786                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1787                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1788                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
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 */
1795              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1796  #pragma omp for nowait  #pragma omp for nowait
1797                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1798                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1799                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
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 */
1806              /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */          } /* end of parallel section */
         } // end of parallel section  
1807      } else {      } else {
1808          out.requireWrite();          out.requireWrite();
1809          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1810          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1811  #pragma omp parallel  #pragma omp parallel
1812          {          {
1813              /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */              vector<double> f_00(numComp);
1814                vector<double> f_01(numComp);
1815                vector<double> f_10(numComp);
1816                vector<double> f_11(numComp);
1817              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1818  #pragma omp for nowait  #pragma omp for nowait
1819                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1820                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1821                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1822                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1823                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1824                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_01[i] + c1*f_00[i];
1825                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_01[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_00[i] + c1*f_01[i];
1826                      } /* end of component loop i */                      } /* end of component loop i */
1827                  } /* end of k1 loop */                  } /* end of k1 loop */
1828              } /* end of face 0 */              } /* end of face 0 */
1829              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1830  #pragma omp for nowait  #pragma omp for nowait
1831                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1832                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1833                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1834                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1835                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1836                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = c1*f_10[i] + c0*f_11[i];
1837                          o[INDEX2(i,numComp,1)] = f_10[i]*c0 + f_11[i]*c1;                          o[INDEX2(i,numComp,1)] = c1*f_11[i] + c0*f_10[i];
1838                      } /* end of component loop i */                      } /* end of component loop i */
1839                  } /* end of k1 loop */                  } /* end of k1 loop */
1840              } /* end of face 1 */              } /* end of face 1 */
1841              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1842  #pragma omp for nowait  #pragma omp for nowait
1843                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1844                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1845                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1846                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1847                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1848                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_10[i] + c1*f_00[i];
1849                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_10[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_00[i] + c1*f_10[i];
1850                      } /* end of component loop i */                      } /* end of component loop i */
1851                  } /* end of k0 loop */                  } /* end of k0 loop */
1852              } /* end of face 2 */              } /* end of face 2 */
1853              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1854  #pragma omp for nowait  #pragma omp for nowait
1855                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1856                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1857                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1858                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1859                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1860                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_11[i] + c1*f_01[i];
1861                          o[INDEX2(i,numComp,1)] = f_01[i]*c0 + f_11[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_01[i] + c1*f_11[i];
1862                      } /* end of component loop i */                      } /* end of component loop i */
1863                  } /* end of k0 loop */                  } /* end of k0 loop */
1864              } /* end of face 3 */              } /* end of face 3 */
1865              /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */          } /* end of parallel section */
         } // end of parallel section  
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 escript::Data& d, const escript::Data& y) const  
 {  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */  
     const double w0 = -0.1555021169820365539*h1/h0;  
     const double w1 = 0.041666666666666666667;  
     const double w10 = -0.041666666666666666667*h0/h1;  
     const double w11 = 0.1555021169820365539*h1/h0;  
     const double w12 = 0.1555021169820365539*h0/h1;  
     const double w13 = 0.01116454968463011277*h0/h1;  
     const double w14 = 0.01116454968463011277*h1/h0;  
     const double w15 = 0.041666666666666666667*h1/h0;  
     const double w16 = -0.01116454968463011277*h0/h1;  
     const double w17 = -0.1555021169820365539*h0/h1;  
     const double w18 = -0.33333333333333333333*h1/h0;  
     const double w19 = 0.25000000000000000000;  
     const double w2 = -0.15550211698203655390;  
     const double w20 = -0.25000000000000000000;  
     const double w21 = 0.16666666666666666667*h0/h1;  
     const double w22 = -0.16666666666666666667*h1/h0;  
     const double w23 = -0.16666666666666666667*h0/h1;  
     const double w24 = 0.33333333333333333333*h1/h0;  
     const double w25 = 0.33333333333333333333*h0/h1;  
     const double w26 = 0.16666666666666666667*h1/h0;  
     const double w27 = -0.33333333333333333333*h0/h1;  
     const double w28 = -0.032861463941450536761*h1;  
     const double w29 = -0.032861463941450536761*h0;  
     const double w3 = 0.041666666666666666667*h0/h1;  
     const double w30 = -0.12264065304058601714*h1;  
     const double w31 = -0.0023593469594139828636*h1;  
     const double w32 = -0.008805202725216129906*h0;  
     const double w33 = -0.008805202725216129906*h1;  
     const double w34 = 0.032861463941450536761*h1;  
     const double w35 = 0.008805202725216129906*h1;  
     const double w36 = 0.008805202725216129906*h0;  
     const double w37 = 0.0023593469594139828636*h1;  
     const double w38 = 0.12264065304058601714*h1;  
     const double w39 = 0.032861463941450536761*h0;  
     const double w4 = 0.15550211698203655390;  
     const double w40 = -0.12264065304058601714*h0;  
     const double w41 = -0.0023593469594139828636*h0;  
     const double w42 = 0.0023593469594139828636*h0;  
     const double w43 = 0.12264065304058601714*h0;  
     const double w44 = -0.16666666666666666667*h1;  
     const double w45 = -0.083333333333333333333*h0;  
     const double w46 = 0.083333333333333333333*h1;  
     const double w47 = 0.16666666666666666667*h1;  
     const double w48 = 0.083333333333333333333*h0;  
     const double w49 = -0.16666666666666666667*h0;  
     const double w5 = -0.041666666666666666667;  
     const double w50 = 0.16666666666666666667*h0;  
     const double w51 = -0.083333333333333333333*h1;  
     const double w52 = 0.025917019497006092316*h0*h1;  
     const double w53 = 0.0018607582807716854616*h0*h1;  
     const double w54 = 0.0069444444444444444444*h0*h1;  
     const double w55 = 0.09672363354357992482*h0*h1;  
     const double w56 = 0.00049858867864229740201*h0*h1;  
     const double w57 = 0.055555555555555555556*h0*h1;  
     const double w58 = 0.027777777777777777778*h0*h1;  
     const double w59 = 0.11111111111111111111*h0*h1;  
     const double w6 = -0.01116454968463011277*h1/h0;  
     const double w60 = -0.19716878364870322056*h1;  
     const double w61 = -0.19716878364870322056*h0;  
     const double w62 = -0.052831216351296779436*h0;  
     const double w63 = -0.052831216351296779436*h1;  
     const double w64 = 0.19716878364870322056*h1;  
     const double w65 = 0.052831216351296779436*h1;  
     const double w66 = 0.19716878364870322056*h0;  
     const double w67 = 0.052831216351296779436*h0;  
     const double w68 = -0.5*h1;  
     const double w69 = -0.5*h0;  
     const double w7 = 0.011164549684630112770;  
     const double w70 = 0.5*h1;  
     const double w71 = 0.5*h0;  
     const double w72 = 0.1555021169820365539*h0*h1;  
     const double w73 = 0.041666666666666666667*h0*h1;  
     const double w74 = 0.01116454968463011277*h0*h1;  
     const double w75 = 0.25*h0*h1;  
     const double w8 = -0.011164549684630112770;  
     const double w9 = -0.041666666666666666667*h1/h0;  
     /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */  
1869    
1870      rhs.requireWrite();  
1871  #pragma omp parallel  namespace
1872    {
1873        // Calculates a guassian blur colvolution matrix for 2D
1874        // See wiki article on the subject    
1875        double* get2DGauss(unsigned radius, double sigma)
1876      {      {
1877          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          double* arr=new double[(radius*2+1)*(radius*2+1)];
1878  #pragma omp for          double common=M_1_PI*0.5*1/(sigma*sigma);
1879              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {          double total=0;
1880                  for (index_t k0=0; k0<m_NE0; ++k0)  {          int r=static_cast<int>(radius);
1881                      bool add_EM_S=false;          for (int y=-r;y<=r;++y)
1882                      bool add_EM_F=false;          {
1883                      vector<double> EM_S(4*4, 0);              for (int x=-r;x<=r;++x)
1884                      vector<double> EM_F(4, 0);              {        
1885                      const index_t e = k0 + m_NE0*k1;                  arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));
1886                      /*** GENERATOR SNIP_PDE_SINGLE TOP */                  total+=arr[(x+r)+(y+r)*(r*2+1)];
1887                      ///////////////              }
1888                      // process A //          }
1889                      ///////////////          double invtotal=1/total;
1890                      if (!A.isEmpty()) {          for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)
1891                          add_EM_S=true;          {
1892                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);              arr[p]*=invtotal;
1893                          if (A.actsExpanded()) {          }
1894                              const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];          return arr;
1895                              const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];      }
1896                              const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];      
1897                              const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];      // applies conv to source to get a point.
1898                              const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];      // (xp, yp) are the coords in the source matrix not the destination matrix
1899                              const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)
1900                              const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];      {
1901                              const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];          size_t bx=xp-radius, by=yp-radius;
1902                              const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];          size_t sbase=bx+by*width;
1903                              const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];          double result=0;
1904                              const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];          for (int y=0;y<2*radius+1;++y)
1905                              const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];          {        
1906                              const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];              for (int x=0;x<2*radius+1;++x)
1907                              const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];              {
1908                              const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];                  result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];
1909                              const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];              }
1910                              const register double tmp4_0 = A_10_1 + A_10_2;          }
1911                              const register double tmp12_0 = A_11_0 + A_11_2;          return result;      
1912                              const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;      }
1913                              const register double tmp10_0 = A_01_3 + A_10_3;  }
                             const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;  
                             const register double tmp0_0 = A_01_0 + A_01_3;  
                             const register double tmp13_0 = A_01_2 + A_10_1;  
                             const register double tmp3_0 = A_00_2 + A_00_3;  
                             const register double tmp11_0 = A_11_1 + A_11_3;  
                             const register double tmp18_0 = A_01_1 + A_10_1;  
                             const register double tmp1_0 = A_00_0 + A_00_1;  
                             const register double tmp15_0 = A_01_1 + A_10_2;  
                             const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;  
                             const register double tmp16_0 = A_10_0 + A_10_3;  
                             const register double tmp6_0 = A_01_3 + A_10_0;  
                             const register double tmp17_0 = A_01_1 + A_01_2;  
                             const register double tmp9_0 = A_01_0 + A_10_0;  
                             const register double tmp7_0 = A_01_0 + A_10_3;  
                             const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;  
                             const register double tmp19_0 = A_01_2 + A_10_2;  
                             const register double tmp14_1 = A_10_0*w8;  
                             const register double tmp23_1 = tmp3_0*w14;  
                             const register double tmp35_1 = A_01_0*w8;  
                             const register double tmp54_1 = tmp13_0*w8;  
                             const register double tmp20_1 = tmp9_0*w4;  
                             const register double tmp25_1 = tmp12_0*w12;  
                             const register double tmp2_1 = A_01_1*w4;  
                             const register double tmp44_1 = tmp7_0*w7;  
                             const register double tmp26_1 = tmp10_0*w4;  
                             const register double tmp52_1 = tmp18_0*w8;  
                             const register double tmp48_1 = A_10_1*w7;  
                             const register double tmp46_1 = A_01_3*w8;  
                             const register double tmp50_1 = A_01_0*w2;  
                             const register double tmp8_1 = tmp4_0*w5;  
                             const register double tmp56_1 = tmp19_0*w8;  
                             const register double tmp9_1 = tmp2_0*w10;  
                             const register double tmp19_1 = A_10_3*w2;  
                             const register double tmp47_1 = A_10_2*w4;  
                             const register double tmp16_1 = tmp3_0*w0;  
                             const register double tmp18_1 = tmp1_0*w6;  
                             const register double tmp31_1 = tmp11_0*w12;  
                             const register double tmp55_1 = tmp15_0*w2;  
                             const register double tmp39_1 = A_10_2*w7;  
                             const register double tmp11_1 = tmp6_0*w7;  
                             const register double tmp40_1 = tmp11_0*w17;  
                             const register double tmp34_1 = tmp15_0*w8;  
                             const register double tmp33_1 = tmp14_0*w5;  
                             const register double tmp24_1 = tmp11_0*w13;  
                             const register double tmp3_1 = tmp1_0*w0;  
                             const register double tmp5_1 = tmp2_0*w3;  
                             const register double tmp43_1 = tmp17_0*w5;  
                             const register double tmp15_1 = A_01_2*w4;  
                             const register double tmp53_1 = tmp19_0*w2;  
                             const register double tmp27_1 = tmp3_0*w11;  
                             const register double tmp32_1 = tmp13_0*w2;  
                             const register double tmp10_1 = tmp5_0*w9;  
                             const register double tmp37_1 = A_10_1*w4;  
                             const register double tmp38_1 = tmp5_0*w15;  
                             const register double tmp17_1 = A_01_1*w7;  
                             const register double tmp12_1 = tmp7_0*w4;  
                             const register double tmp22_1 = tmp10_0*w7;  
                             const register double tmp57_1 = tmp18_0*w2;  
                             const register double tmp28_1 = tmp9_0*w7;  
                             const register double tmp29_1 = tmp1_0*w14;  
                             const register double tmp51_1 = tmp11_0*w16;  
                             const register double tmp42_1 = tmp12_0*w16;  
                             const register double tmp49_1 = tmp12_0*w17;  
                             const register double tmp21_1 = tmp1_0*w11;  
                             const register double tmp1_1 = tmp0_0*w1;  
                             const register double tmp45_1 = tmp6_0*w4;  
                             const register double tmp7_1 = A_10_0*w2;  
                             const register double tmp6_1 = tmp3_0*w6;  
                             const register double tmp13_1 = tmp8_0*w1;  
                             const register double tmp36_1 = tmp16_0*w1;  
                             const register double tmp41_1 = A_01_3*w2;  
                             const register double tmp30_1 = tmp12_0*w13;  
                             const register double tmp4_1 = A_01_2*w7;  
                             const register double tmp0_1 = A_10_3*w8;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
                             EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;  
                         } else { /* constant data */  
                             const register double A_00 = A_p[INDEX2(0,0,2)];  
                             const register double A_01 = A_p[INDEX2(0,1,2)];  
                             const register double A_10 = A_p[INDEX2(1,0,2)];  
                             const register double A_11 = A_p[INDEX2(1,1,2)];  
                             const register double tmp0_0 = A_01 + A_10;  
                             const register double tmp0_1 = A_00*w18;  
                             const register double tmp10_1 = A_01*w20;  
                             const register double tmp12_1 = A_00*w26;  
                             const register double tmp4_1 = A_00*w22;  
                             const register double tmp8_1 = A_00*w24;  
                             const register double tmp13_1 = A_10*w19;  
                             const register double tmp9_1 = tmp0_0*w20;  
                             const register double tmp3_1 = A_11*w21;  
                             const register double tmp11_1 = A_11*w27;  
                             const register double tmp1_1 = A_01*w19;  
                             const register double tmp6_1 = A_11*w23;  
                             const register double tmp7_1 = A_11*w25;  
                             const register double tmp2_1 = A_10*w20;  
                             const register double tmp5_1 = tmp0_0*w19;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                         }  
                     }  
                     ///////////////  
                     // process B //  
                     ///////////////  
                     if (!B.isEmpty()) {  
                         add_EM_S=true;  
                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);  
                         if (B.actsExpanded()) {  
                             const register double B_0_0 = B_p[INDEX2(0,0,2)];  
                             const register double B_1_0 = B_p[INDEX2(1,0,2)];  
                             const register double B_0_1 = B_p[INDEX2(0,1,2)];  
                             const register double B_1_1 = B_p[INDEX2(1,1,2)];  
                             const register double B_0_2 = B_p[INDEX2(0,2,2)];  
                             const register double B_1_2 = B_p[INDEX2(1,2,2)];  
                             const register double B_0_3 = B_p[INDEX2(0,3,2)];  
                             const register double B_1_3 = B_p[INDEX2(1,3,2)];  
                             const register double tmp3_0 = B_0_0 + B_0_2;  
                             const register double tmp1_0 = B_1_2 + B_1_3;  
                             const register double tmp2_0 = B_0_1 + B_0_3;  
                             const register double tmp0_0 = B_1_0 + B_1_1;  
                             const register double tmp63_1 = B_1_1*w42;  
                             const register double tmp79_1 = B_1_1*w40;  
                             const register double tmp37_1 = tmp3_0*w35;  
                             const register double tmp8_1 = tmp0_0*w32;  
                             const register double tmp71_1 = B_0_1*w34;  
                             const register double tmp19_1 = B_0_3*w31;  
                             const register double tmp15_1 = B_0_3*w34;  
                             const register double tmp9_1 = tmp3_0*w34;  
                             const register double tmp35_1 = B_1_0*w36;  
                             const register double tmp66_1 = B_0_3*w28;  
                             const register double tmp28_1 = B_1_0*w42;  
                             const register double tmp22_1 = B_1_0*w40;  
                             const register double tmp16_1 = B_1_2*w29;  
                             const register double tmp6_1 = tmp2_0*w35;  
                             const register double tmp55_1 = B_1_3*w40;  
                             const register double tmp50_1 = B_1_3*w42;  
                             const register double tmp7_1 = tmp1_0*w29;  
                             const register double tmp1_1 = tmp1_0*w32;  
                             const register double tmp57_1 = B_0_3*w30;  
                             const register double tmp18_1 = B_1_1*w32;  
                             const register double tmp53_1 = B_1_0*w41;  
                             const register double tmp61_1 = B_1_3*w36;  
                             const register double tmp27_1 = B_0_3*w38;  
                             const register double tmp64_1 = B_0_2*w30;  
                             const register double tmp76_1 = B_0_1*w38;  
                             const register double tmp39_1 = tmp2_0*w34;  
                             const register double tmp62_1 = B_0_1*w31;  
                             const register double tmp56_1 = B_0_0*w31;  
                             const register double tmp49_1 = B_1_1*w36;  
                             const register double tmp2_1 = B_0_2*w31;  
                             const register double tmp23_1 = B_0_2*w33;  
                             const register double tmp38_1 = B_1_1*w43;  
                             const register double tmp74_1 = B_1_2*w41;  
                             const register double tmp43_1 = B_1_1*w41;  
                             const register double tmp58_1 = B_0_2*w28;  
                             const register double tmp67_1 = B_0_0*w33;  
                             const register double tmp33_1 = tmp0_0*w39;  
                             const register double tmp4_1 = B_0_0*w28;  
                             const register double tmp20_1 = B_0_0*w30;  
                             const register double tmp13_1 = B_0_2*w38;  
                             const register double tmp65_1 = B_1_2*w43;  
                             const register double tmp0_1 = tmp0_0*w29;  
                             const register double tmp41_1 = tmp3_0*w33;  
                             const register double tmp73_1 = B_0_2*w37;  
                             const register double tmp69_1 = B_0_0*w38;  
                             const register double tmp48_1 = B_1_2*w39;  
                             const register double tmp59_1 = B_0_1*w33;  
                             const register double tmp17_1 = B_1_3*w41;  
                             const register double tmp5_1 = B_0_3*w33;  
                             const register double tmp3_1 = B_0_1*w30;  
                             const register double tmp21_1 = B_0_1*w28;  
                             const register double tmp42_1 = B_1_0*w29;  
                             const register double tmp54_1 = B_1_2*w32;  
                             const register double tmp60_1 = B_1_0*w39;  
                             const register double tmp32_1 = tmp1_0*w36;  
                             const register double tmp10_1 = B_0_1*w37;  
                             const register double tmp14_1 = B_0_0*w35;  
                             const register double tmp29_1 = B_0_1*w35;  
                             const register double tmp26_1 = B_1_2*w36;  
                             const register double tmp30_1 = B_1_3*w43;  
                             const register double tmp70_1 = B_0_2*w35;  
                             const register double tmp34_1 = B_1_3*w39;  
                             const register double tmp51_1 = B_1_0*w43;  
                             const register double tmp31_1 = B_0_2*w34;  
                             const register double tmp45_1 = tmp3_0*w28;  
                             const register double tmp11_1 = tmp1_0*w39;  
                             const register double tmp52_1 = B_1_1*w29;  
                             const register double tmp44_1 = B_1_3*w32;  
                             const register double tmp25_1 = B_1_1*w39;  
                             const register double tmp47_1 = tmp2_0*w33;  
                             const register double tmp72_1 = B_1_3*w29;  
                             const register double tmp40_1 = tmp2_0*w28;  
                             const register double tmp46_1 = B_1_2*w40;  
                             const register double tmp36_1 = B_1_2*w42;  
                             const register double tmp24_1 = B_0_0*w37;  
                             const register double tmp77_1 = B_0_3*w35;  
                             const register double tmp68_1 = B_0_3*w37;  
                             const register double tmp78_1 = B_0_0*w34;  
                             const register double tmp12_1 = tmp0_0*w36;  
                             const register double tmp75_1 = B_1_0*w32;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                             EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                         } else { /* constant data */  
                             const register double B_0 = B_p[0];  
                             const register double B_1 = B_p[1];  
                             const register double tmp6_1 = B_1*w50;  
                             const register double tmp1_1 = B_1*w45;  
                             const register double tmp5_1 = B_1*w49;  
                             const register double tmp4_1 = B_1*w48;  
                             const register double tmp0_1 = B_0*w44;  
                             const register double tmp2_1 = B_0*w46;  
                             const register double tmp7_1 = B_0*w51;  
                             const register double tmp3_1 = B_0*w47;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;  
                         }  
                     }  
                     ///////////////  
                     // process C //  
                     ///////////////  
                     if (!C.isEmpty()) {  
                         add_EM_S=true;  
                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);  
                         if (C.actsExpanded()) {  
                             const register double C_0_0 = C_p[INDEX2(0,0,2)];  
                             const register double C_1_0 = C_p[INDEX2(1,0,2)];  
                             const register double C_0_1 = C_p[INDEX2(0,1,2)];  
                             const register double C_1_1 = C_p[INDEX2(1,1,2)];  
                             const register double C_0_2 = C_p[INDEX2(0,2,2)];  
                             const register double C_1_2 = C_p[INDEX2(1,2,2)];  
                             const register double C_0_3 = C_p[INDEX2(0,3,2)];  
                             const register double C_1_3 = C_p[INDEX2(1,3,2)];  
                             const register double tmp2_0 = C_0_1 + C_0_3;  
                             const register double tmp1_0 = C_1_2 + C_1_3;  
                             const register double tmp3_0 = C_0_0 + C_0_2;  
                             const register double tmp0_0 = C_1_0 + C_1_1;  
                             const register double tmp64_1 = C_0_2*w30;  
                             const register double tmp14_1 = C_0_2*w28;  
                             const register double tmp19_1 = C_0_3*w31;  
                             const register double tmp22_1 = C_1_0*w40;  
                             const register double tmp37_1 = tmp3_0*w35;  
                             const register double tmp29_1 = C_0_1*w35;  
                             const register double tmp73_1 = C_0_2*w37;  
                             const register double tmp74_1 = C_1_2*w41;  
                             const register double tmp52_1 = C_1_3*w39;  
                             const register double tmp25_1 = C_1_1*w39;  
                             const register double tmp62_1 = C_0_1*w31;  
                             const register double tmp79_1 = C_1_1*w40;  
                             const register double tmp43_1 = C_1_1*w36;  
                             const register double tmp27_1 = C_0_3*w38;  
                             const register double tmp28_1 = C_1_0*w42;  
                             const register double tmp63_1 = C_1_1*w42;  
                             const register double tmp59_1 = C_0_3*w34;  
                             const register double tmp72_1 = C_1_3*w29;  
                             const register double tmp40_1 = tmp2_0*w35;  
                             const register double tmp13_1 = C_0_3*w30;  
                             const register double tmp51_1 = C_1_2*w40;  
                             const register double tmp54_1 = C_1_2*w42;  
                             const register double tmp12_1 = C_0_0*w31;  
                             const register double tmp2_1 = tmp1_0*w32;  
                             const register double tmp68_1 = C_0_2*w31;  
                             const register double tmp75_1 = C_1_0*w32;  
                             const register double tmp49_1 = C_1_1*w41;  
                             const register double tmp4_1 = C_0_2*w35;  
                             const register double tmp66_1 = C_0_3*w28;  
                             const register double tmp56_1 = C_0_1*w37;  
                             const register double tmp5_1 = C_0_1*w34;  
                             const register double tmp38_1 = tmp2_0*w34;  
                             const register double tmp76_1 = C_0_1*w38;  
                             const register double tmp21_1 = C_0_1*w28;  
                             const register double tmp69_1 = C_0_1*w30;  
                             const register double tmp53_1 = C_1_0*w36;  
                             const register double tmp42_1 = C_1_2*w39;  
                             const register double tmp32_1 = tmp1_0*w29;  
                             const register double tmp45_1 = C_1_0*w43;  
                             const register double tmp33_1 = tmp0_0*w32;  
                             const register double tmp35_1 = C_1_0*w41;  
                             const register double tmp26_1 = C_1_2*w36;  
                             const register double tmp67_1 = C_0_0*w33;  
                             const register double tmp31_1 = C_0_2*w34;  
                             const register double tmp20_1 = C_0_0*w30;  
                             const register double tmp70_1 = C_0_0*w28;  
                             const register double tmp8_1 = tmp0_0*w39;  
                             const register double tmp30_1 = C_1_3*w43;  
                             const register double tmp0_1 = tmp0_0*w29;  
                             const register double tmp17_1 = C_1_3*w41;  
                             const register double tmp58_1 = C_0_0*w35;  
                             const register double tmp9_1 = tmp3_0*w33;  
                             const register double tmp61_1 = C_1_3*w36;  
                             const register double tmp41_1 = tmp3_0*w34;  
                             const register double tmp50_1 = C_1_3*w32;  
                             const register double tmp18_1 = C_1_1*w32;  
                             const register double tmp6_1 = tmp1_0*w36;  
                             const register double tmp3_1 = C_0_0*w38;  
                             const register double tmp34_1 = C_1_1*w29;  
                             const register double tmp77_1 = C_0_3*w35;  
                             const register double tmp65_1 = C_1_2*w43;  
                             const register double tmp71_1 = C_0_3*w33;  
                             const register double tmp55_1 = C_1_1*w43;  
                             const register double tmp46_1 = tmp3_0*w28;  
                             const register double tmp24_1 = C_0_0*w37;  
                             const register double tmp10_1 = tmp1_0*w39;  
                             const register double tmp48_1 = C_1_0*w29;  
                             const register double tmp15_1 = C_0_1*w33;  
                             const register double tmp36_1 = C_1_2*w32;  
                             const register double tmp60_1 = C_1_0*w39;  
                             const register double tmp47_1 = tmp2_0*w33;  
                             const register double tmp16_1 = C_1_2*w29;  
                             const register double tmp1_1 = C_0_3*w37;  
                             const register double tmp7_1 = tmp2_0*w28;  
                             const register double tmp39_1 = C_1_3*w40;  
                             const register double tmp44_1 = C_1_3*w42;  
                             const register double tmp57_1 = C_0_2*w38;  
                             const register double tmp78_1 = C_0_0*w34;  
                             const register double tmp11_1 = tmp0_0*w36;  
                             const register double tmp23_1 = C_0_2*w33;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                             EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                         } else { /* constant data */  
                             const register double C_0 = C_p[0];  
                             const register double C_1 = C_p[1];  
                             const register double tmp1_1 = C_1*w45;  
                             const register double tmp3_1 = C_0*w51;  
                             const register double tmp4_1 = C_0*w44;  
                             const register double tmp7_1 = C_0*w46;  
                             const register double tmp5_1 = C_1*w49;  
                             const register double tmp2_1 = C_1*w48;  
                             const register double tmp0_1 = C_0*w47;  
                             const register double tmp6_1 = C_1*w50;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;  
                             EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;  
                         }  
                     }  
                     ///////////////  
                     // process D //  
                     ///////////////  
                     if (!D.isEmpty()) {  
                         add_EM_S=true;  
                         const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);  
                         if (D.actsExpanded()) {  
                             const register double D_0 = D_p[0];  
                             const register double D_1 = D_p[1];  
                             const register double D_2 = D_p[2];  
                             const register double D_3 = D_p[3];  
                             const register double tmp4_0 = D_1 + D_3;  
                             const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;  
                             const register double tmp5_0 = D_0 + D_2;  
                             const register double tmp0_0 = D_0 + D_1;  
                             const register double tmp6_0 = D_0 + D_3;  
                             const register double tmp1_0 = D_2 + D_3;  
                             const register double tmp3_0 = D_1 + D_2;  
                             const register double tmp16_1 = D_1*w56;  
                             const register double tmp14_1 = tmp6_0*w54;  
                             const register double tmp8_1 = D_3*w55;  
                             const register double tmp2_1 = tmp2_0*w54;  
                             const register double tmp12_1 = tmp5_0*w52;  
                             const register double tmp4_1 = tmp0_0*w53;  
                             const register double tmp3_1 = tmp1_0*w52;  
                             const register double tmp13_1 = tmp4_0*w53;  
                             const register double tmp10_1 = tmp4_0*w52;  
                             const register double tmp1_1 = tmp1_0*w53;  
                             const register double tmp7_1 = D_3*w56;  
                             const register double tmp0_1 = tmp0_0*w52;  
                             const register double tmp11_1 = tmp5_0*w53;  
                             const register double tmp9_1 = D_0*w56;  
                             const register double tmp5_1 = tmp3_0*w54;  
                             const register double tmp18_1 = D_2*w56;  
                             const register double tmp17_1 = D_1*w55;  
                             const register double tmp6_1 = D_0*w55;  
                             const register double tmp15_1 = D_2*w55;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp2_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp2_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;  
                         } else { /* constant data */  
                             const register double D_0 = D_p[0];  
                             const register double tmp2_1 = D_0*w59;  
                             const register double tmp1_1 = D_0*w58;  
                             const register double tmp0_1 = D_0*w57;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1;  
                             EM_S[INDEX2(0,0,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp1_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp2_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp1_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp2_1;  
                         }  
                     }  
                     ///////////////  
                     // process X //  
                     ///////////////  
                     if (!X.isEmpty()) {  
                         add_EM_F=true;  
                         const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);  
                         if (X.actsExpanded()) {  
                             const register double X_0_0 = X_p[INDEX2(0,0,2)];  
                             const register double X_1_0 = X_p[INDEX2(1,0,2)];  
                             const register double X_0_1 = X_p[INDEX2(0,1,2)];  
                             const register double X_1_1 = X_p[INDEX2(1,1,2)];  
                             const register double X_0_2 = X_p[INDEX2(0,2,2)];  
                             const register double X_1_2 = X_p[INDEX2(1,2,2)];  
                             const register double X_0_3 = X_p[INDEX2(0,3,2)];  
                             const register double X_1_3 = X_p[INDEX2(1,3,2)];  
                             const register double tmp1_0 = X_1_1 + X_1_3;  
                             const register double tmp3_0 = X_0_0 + X_0_1;  
                             const register double tmp2_0 = X_1_0 + X_1_2;  
                             const register double tmp0_0 = X_0_2 + X_0_3;  
                             const register double tmp8_1 = tmp2_0*w66;  
                             const register double tmp5_1 = tmp3_0*w64;  
                             const register double tmp14_1 = tmp0_0*w64;  
                             const register double tmp3_1 = tmp3_0*w60;  
                             const register double tmp9_1 = tmp3_0*w63;  
                             const register double tmp13_1 = tmp3_0*w65;  
                             const register double tmp12_1 = tmp1_0*w66;  
                             const register double tmp10_1 = tmp0_0*w60;  
                             const register double tmp2_1 = tmp2_0*w61;  
                             const register double tmp6_1 = tmp2_0*w62;  
                             const register double tmp4_1 = tmp0_0*w65;  
                             const register double tmp11_1 = tmp1_0*w67;  
                             const register double tmp1_1 = tmp1_0*w62;  
                             const register double tmp7_1 = tmp1_0*w61;  
                             const register double tmp0_1 = tmp0_0*w63;  
                             const register double tmp15_1 = tmp2_0*w67;  
                             EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;  
                             EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;  
                             EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                         } else { /* constant data */  
                             const register double X_0 = X_p[0];  
                             const register double X_1 = X_p[1];  
                             const register double tmp3_1 = X_1*w71;  
                             const register double tmp2_1 = X_0*w70;  
                             const register double tmp1_1 = X_0*w68;  
                             const register double tmp0_1 = X_1*w69;  
                             EM_F[0]+=tmp0_1 + tmp1_1;  
                             EM_F[1]+=tmp0_1 + tmp2_1;  
                             EM_F[2]+=tmp1_1 + tmp3_1;  
                             EM_F[3]+=tmp2_1 + tmp3_1;  
                         }  
                     }  
                     ///////////////  
                     // process Y //  
                     ///////////////  
                     if (!Y.isEmpty()) {  
                         add_EM_F=true;  
                         const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);  
                         if (Y.actsExpanded()) {  
                             const register double Y_0 = Y_p[0];  
                             const register double Y_1 = Y_p[1];  
                             const register double Y_2 = Y_p[2];  
                             const register double Y_3 = Y_p[3];  
                             const register double tmp0_0 = Y_1 + Y_2;  
                             const register double tmp1_0 = Y_0 + Y_3;  
                             const register double tmp9_1 = Y_0*w74;  
                             const register double tmp4_1 = tmp1_0*w73;  
                             const register double tmp5_1 = Y_2*w74;  
                             const register double tmp7_1 = Y_1*w74;  
                             const register double tmp6_1 = Y_2*w72;  
                             const register double tmp2_1 = Y_3*w74;  
                             const register double tmp8_1 = Y_3*w72;  
                             const register double tmp3_1 = Y_1*w72;  
                             const register double tmp0_1 = Y_0*w72;  
                             const register double tmp1_1 = tmp0_0*w73;  
                             EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;  
                             EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;  
                             EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;  
                         } else { /* constant data */  
                             const register double Y_0 = Y_p[0];  
                             const register double tmp0_1 = Y_0*w75;  
                             EM_F[0]+=tmp0_1;  
                             EM_F[1]+=tmp0_1;  
                             EM_F[2]+=tmp0_1;  
                             EM_F[3]+=tmp0_1;  
                         }  
                     }  
                     /* GENERATOR SNIP_PDE_SINGLE BOTTOM */  
1914    
                     // add to matrix (if add_EM_S) and RHS (if add_EM_F)  
                     const index_t firstNode=m_N0*k1+k0;  
                     IndexVector rowIndex;  
                     rowIndex.push_back(m_dofMap[firstNode]);  
                     rowIndex.push_back(m_dofMap[firstNode+1]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                     if (add_EM_F) {  
                         //cout << "-- AddtoRHS -- " << endl;  
                         double *F_p=rhs.getSampleDataRW(0);  
                         for (index_t i=0; i<rowIndex.size(); i++) {  
                             if (rowIndex[i]<getNumDOF()) {  
                                 F_p[rowIndex[i]]+=EM_F[i];  
                                 //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;  
                             }  
                         }  
                         //cout << "---"<<endl;  
                     }  
                     if (add_EM_S) {  
                         //cout << "-- AddtoSystem -- " << endl;  
                         addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);  
                     }  
1915    
1916                  } // end k0 loop  /* This is a wrapper for filtered (and non-filtered) randoms
1917              } // end k1 loop   * For detailed doco see randomFillWorker
1918          } // end of colouring  */
1919      } // end of parallel region  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            throw RipleyException("Ripley only supports filters for scalar data.");
1927        }
1928        escript::Data res=randomFillWorker(shape, seed, filter);
1929        if (res.getFunctionSpace()!=what)
1930        {
1931            escript::Data r=escript::Data(res, what);
1932            return r;
1933        }
1934        return res;
1935  }  }
1936    
 //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 escript::Data& d, const escript::Data& y) const  
 {  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     /*** GENERATOR SNIP_PDE_SINGLE_REDUCED_PRE TOP */  
     const double w0 = -0.25*h1/h0;  
     const double w1 = 0.25;  
     const double w2 = -0.25;  
     const double w3 = 0.25*h0/h1;  
     const double w4 = -0.25*h0/h1;  
     const double w5 = 0.25*h1/h0;  
     const double w6 = -0.125*h1;  
     const double w7 = -0.125*h0;  
     const double w8 = 0.125*h1;  
     const double w9 = 0.125*h0;  
     const double w10 = 0.0625*h0*h1;  
     const double w11 = -0.5*h1;  
     const double w12 = -0.5*h0;  
     const double w13 = 0.5*h1;  
     const double w14 = 0.5*h0;  
     const double w15 = 0.25*h0*h1;  
     /* GENERATOR SNIP_PDE_SINGLE_REDUCED_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_NE1; k1+=2) {  
                 for (index_t k0=0; k0<m_NE0; ++k0)  {  
                     bool add_EM_S=false;  
                     bool add_EM_F=false;  
                     vector<double> EM_S(4*4, 0);  
                     vector<double> EM_F(4, 0);  
                     const index_t e = k0 + m_NE0*k1;  
                     /*** GENERATOR SNIP_PDE_SINGLE_REDUCED TOP */  
                     ///////////////  
                     // process A //  
                     ///////////////  
                     if (!A.isEmpty()) {  
                         add_EM_S=true;  
                         const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);  
                         const register double A_00 = A_p[INDEX2(0,0,2)];  
                         const register double A_01 = A_p[INDEX2(0,1,2)];  
                         const register double A_10 = A_p[INDEX2(1,0,2)];  
                         const register double A_11 = A_p[INDEX2(1,1,2)];  
                         const register double tmp0_0 = A_01 + A_10;  
                         const register double tmp2_1 = A_01*w1;  
                         const register double tmp7_1 = tmp0_0*w2;  
                         const register double tmp3_1 = A_10*w2;  
                         const register double tmp6_1 = A_00*w5;  
                         const register double tmp1_1 = A_00*w0;  
                         const register double tmp4_1 = tmp0_0*w1;  
                         const register double tmp9_1 = A_01*w2;  
                         const register double tmp5_1 = A_11*w4;  
                         const register double tmp8_1 = A_10*w1;  
                         const register double tmp0_1 = A_11*w3;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;  
                         EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;  
                         EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;  
                         EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;  
                         EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;  
                         EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
                         EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;  
                         EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;  
                         EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;  
                     }  
                     ///////////////  
                     // process B //  
                     ///////////////  
                     if (!B.isEmpty()) {  
                         add_EM_S=true;  
                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);  
                         const register double B_0 = B_p[0];  
                         const register double B_1 = B_p[1];  
                         const register double tmp2_1 = B_0*w8;  
                         const register double tmp0_1 = B_0*w6;  
                         const register double tmp3_1 = B_1*w9;  
                         const register double tmp1_1 = B_1*w7;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;  
                     }  
                     ///////////////  
                     // process C //  
                     ///////////////  
                     if (!C.isEmpty()) {  
                         add_EM_S=true;  
                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);  
                         const register double C_0 = C_p[0];  
                         const register double C_1 = C_p[1];  
                         const register double tmp1_1 = C_1*w7;  
                         const register double tmp0_1 = C_0*w8;  
                         const register double tmp3_1 = C_0*w6;  
                         const register double tmp2_1 = C_1*w9;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;  
                         EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;  
                     }  
                     ///////////////  
                     // process D //  
                     ///////////////  
                     if (!D.isEmpty()) {  
                         add_EM_S=true;  
                         const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);  
                         const register double D_0 = D_p[0];  
                         const register double tmp0_1 = D_0*w10;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp0_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp0_1;  
                         EM_S[INDEX2(0,0,4)]+=tmp0_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1;  
                         EM_S[INDEX2(3,0,4)]+=tmp0_1;  
                         EM_S[INDEX2(3,1,4)]+=tmp0_1;  
                         EM_S[INDEX2(2,1,4)]+=tmp0_1;  
                         EM_S[INDEX2(0,2,4)]+=tmp0_1;  
                         EM_S[INDEX2(2,0,4)]+=tmp0_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp0_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1;  
                         EM_S[INDEX2(2,2,4)]+=tmp0_1;  
                         EM_S[INDEX2(1,0,4)]+=tmp0_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp0_1;  
                         EM_S[INDEX2(1,1,4)]+=tmp0_1;  
                     }  
                     ///////////////  
                     // process X //  
                     ///////////////  
                     if (!X.isEmpty()) {  
                         add_EM_F=true;  
                         const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);  
                         const register double X_0 = X_p[0];  
                         const register double X_1 = X_p[1];  
                         const register double tmp0_1 = X_0*w11;  
                         const register double tmp2_1 = X_0*w13;  
                         const register double tmp1_1 = X_1*w12;  
                         const register double tmp3_1 = X_1*w14;  
                         EM_F[0]+=tmp0_1 + tmp1_1;  
                         EM_F[1]+=tmp1_1 + tmp2_1;  
                         EM_F[2]+=tmp0_1 + tmp3_1;  
                         EM_F[3]+=tmp2_1 + tmp3_1;  
                     }  
                     ///////////////  
                     // process Y //  
                     ///////////////  
                     if (!Y.isEmpty()) {  
                         add_EM_F=true;  
                         const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);  
                         const register double Y_0 = Y_p[0];  
                         const register double tmp0_1 = Y_0*w15;  
                         EM_F[0]+=tmp0_1;  
                         EM_F[1]+=tmp0_1;  
                         EM_F[2]+=tmp0_1;  
                         EM_F[3]+=tmp0_1;  
                     }  
                     /* GENERATOR SNIP_PDE_SINGLE_REDUCED 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_N0*k1+k0;  dimensions = ext[0] * ext[1].
1946                      IndexVector rowIndex;  Where ext[i]= internal[i]+2*radius.
                     rowIndex.push_back(m_dofMap[firstNode]);  
                     rowIndex.push_back(m_dofMap[firstNode+1]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                     if (add_EM_F) {  
                         //cout << "-- AddtoRHS -- " << endl;  
                         double *F_p=rhs.getSampleDataRW(0);  
                         for (index_t i=0; i<rowIndex.size(); i++) {  
                             if (rowIndex[i]<getNumDOF()) {  
                                 F_p[rowIndex[i]]+=EM_F[i];  
                                 //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;  
                             }  
                         }  
                         //cout << "---"<<endl;  
                     }  
                     if (add_EM_S) {  
                         //cout << "-- AddtoSystem -- " << endl;  
                         addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);  
                     }  
1947    
1948                  } // end k0 loop  Now for MPI there is overlap to deal with. We need to share both the overlapping
1949              } // end k1 loop  values themselves but also the external region.
         } // end of colouring  
     } // end of parallel region  
 }  
1950    
1951  //protected  In a hypothetical 1-D case:
 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 escript::Data& d, const escript::Data& y) const  
 {  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     dim_t numEq, numComp;  
     if (!mat)  
         numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());  
     else {  
         numEq=mat->logical_row_block_size;  
         numComp=mat->logical_col_block_size;  
     }  
   
     /*** GENERATOR SNIP_PDE_SYSTEM_PRE TOP */  
     const double w0 = -0.1555021169820365539*h1/h0;  
     const double w1 = 0.041666666666666666667;  
     const double w10 = -0.041666666666666666667*h0/h1;  
     const double w11 = 0.1555021169820365539*h1/h0;  
     const double w12 = 0.1555021169820365539*h0/h1;  
     const double w13 = 0.01116454968463011277*h0/h1;  
     const double w14 = 0.01116454968463011277*h1/h0;  
     const double w15 = 0.041666666666666666667*h1/h0;  
     const double w16 = -0.01116454968463011277*h0/h1;  
     const double w17 = -0.1555021169820365539*h0/h1;  
     const double w18 = -0.33333333333333333333*h1/h0;  
     const double w19 = 0.25000000000000000000;  
     const double w2 = -0.15550211698203655390;  
     const double w20 = -0.25000000000000000000;  
     const double w21 = 0.16666666666666666667*h0/h1;  
     const double w22 = -0.16666666666666666667*h1/h0;  
     const double w23 = -0.16666666666666666667*h0/h1;  
     const double w24 = 0.33333333333333333333*h1/h0;  
     const double w25 = 0.33333333333333333333*h0/h1;  
     const double w26 = 0.16666666666666666667*h1/h0;  
     const double w27 = -0.33333333333333333333*h0/h1;  
     const double w28 = -0.032861463941450536761*h1;  
     const double w29 = -0.032861463941450536761*h0;  
     const double w3 = 0.041666666666666666667*h0/h1;  
     const double w30 = -0.12264065304058601714*h1;  
     const double w31 = -0.0023593469594139828636*h1;  
     const double w32 = -0.008805202725216129906*h0;  
     const double w33 = -0.008805202725216129906*h1;  
     const double w34 = 0.032861463941450536761*h1;  
     const double w35 = 0.008805202725216129906*h1;  
     const double w36 = 0.008805202725216129906*h0;  
     const double w37 = 0.0023593469594139828636*h1;  
     const double w38 = 0.12264065304058601714*h1;  
     const double w39 = 0.032861463941450536761*h0;  
     const double w4 = 0.15550211698203655390;  
     const double w40 = -0.12264065304058601714*h0;  
     const double w41 = -0.0023593469594139828636*h0;  
     const double w42 = 0.0023593469594139828636*h0;  
     const double w43 = 0.12264065304058601714*h0;  
     const double w44 = -0.16666666666666666667*h1;  
     const double w45 = -0.083333333333333333333*h0;  
     const double w46 = 0.083333333333333333333*h1;  
     const double w47 = 0.16666666666666666667*h1;  
     const double w48 = 0.083333333333333333333*h0;  
     const double w49 = -0.16666666666666666667*h0;  
     const double w5 = -0.041666666666666666667;  
     const double w50 = 0.16666666666666666667*h0;  
     const double w51 = -0.083333333333333333333*h1;  
     const double w52 = 0.025917019497006092316*h0*h1;  
     const double w53 = 0.0018607582807716854616*h0*h1;  
     const double w54 = 0.0069444444444444444444*h0*h1;  
     const double w55 = 0.09672363354357992482*h0*h1;  
     const double w56 = 0.00049858867864229740201*h0*h1;  
     const double w57 = 0.055555555555555555556*h0*h1;  
     const double w58 = 0.027777777777777777778*h0*h1;  
     const double w59 = 0.11111111111111111111*h0*h1;  
     const double w6 = -0.01116454968463011277*h1/h0;  
     const double w60 = -0.19716878364870322056*h1;  
     const double w61 = -0.19716878364870322056*h0;  
     const double w62 = -0.052831216351296779436*h0;  
     const double w63 = -0.052831216351296779436*h1;  
     const double w64 = 0.19716878364870322056*h1;  
     const double w65 = 0.052831216351296779436*h1;  
     const double w66 = 0.19716878364870322056*h0;  
     const double w67 = 0.052831216351296779436*h0;  
     const double w68 = -0.5*h1;  
     const double w69 = -0.5*h0;  
     const double w7 = 0.011164549684630112770;  
     const double w70 = 0.5*h1;  
     const double w71 = 0.5*h0;  
     const double w72 = 0.1555021169820365539*h0*h1;  
     const double w73 = 0.041666666666666666667*h0*h1;  
     const double w74 = 0.01116454968463011277*h0*h1;  
     const double w75 = 0.25*h0*h1;  
     const double w8 = -0.011164549684630112770;  
     const double w9 = -0.041666666666666666667*h1/h0;  
     /* GENERATOR SNIP_PDE_SYSTEM_PRE BOTTOM */  
1952    
1953      rhs.requireWrite();  
1954  #pragma omp parallel  1234567
1955    would be split into two ranks thus:
1956    123(4)  (4)567     [4 being a shared element]
1957    
1958    If the radius is 2.   There will be padding elements on the outside:
1959    
1960    pp123(4)  (4)567pp
1961    
1962    To ensure that 4 can be correctly computed on both ranks, values from the other rank
1963    need to be known.
1964    
1965    pp123(4)56   23(4)567pp
1966    
1967    Now in our case, we wout set all the values 23456 on the left rank and send them to the
1968    right hand rank.
1969    
1970    So the edges _may_ need to be shared at a distance `inset` from all boundaries.
1971    inset=2*radius+1    
1972    This is to ensure that values at distance `radius` from the shared/overlapped element
1973    that ripley has.
1974    
1975    
1976    */
1977    escript::Data Rectangle::randomFillWorker(const escript::DataTypes::ShapeType& shape,
1978           long seed, const boost::python::tuple& filter) const
1979    {
1980        if (m_numDim!=2)
1981      {      {
1982          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          throw RipleyException("Only 2D supported at this time.");
1983  #pragma omp for      }
             for (index_t k1=k1_0; k1<m_NE1; k1+=2) {  
                 for (index_t k0=0; k0<m_NE0; ++k0)  {  
                     bool add_EM_S=false;  
                     bool add_EM_F=false;  
                     vector<double> EM_S(4*4*numEq*numComp, 0);  
                     vector<double> EM_F(4*numEq, 0);  
                     const index_t e = k0 + m_NE0*k1;  
                     /*** 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 register double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];  
                                     const register double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];  
                                     const register double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];  
                                     const register double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];  
                                     const register double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];  
                                     const register double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];  
                                     const register double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];  
                                     const register double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];  
                                     const register double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];  
                                     const register double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];  
                                     const register double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];  
                                     const register double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];  
                                     const register double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];  
                                     const register double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];  
                                     const register double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];  
                                     const register double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];  
                                     const register double tmp4_0 = A_10_1 + A_10_2;  
                                     const register double tmp12_0 = A_11_0 + A_11_2;  
                                     const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;  
                                     const register double tmp10_0 = A_01_3 + A_10_3;  
                                     const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;  
                                     const register double tmp0_0 = A_01_0 + A_01_3;  
                                     const register double tmp13_0 = A_01_2 + A_10_1;  
                                     const register double tmp3_0 = A_00_2 + A_00_3;  
                                     const register double tmp11_0 = A_11_1 + A_11_3;  
                                     const register double tmp18_0 = A_01_1 + A_10_1;  
                                     const register double tmp1_0 = A_00_0 + A_00_1;  
                                     const register double tmp15_0 = A_01_1 + A_10_2;  
                                     const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;  
                                     const register double tmp16_0 = A_10_0 + A_10_3;  
                                     const register double tmp6_0 = A_01_3 + A_10_0;  
                                     const register double tmp17_0 = A_01_1 + A_01_2;  
                                     const register double tmp9_0 = A_01_0 + A_10_0;  
                                     const register double tmp7_0 = A_01_0 + A_10_3;  
                                     const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;  
                                     const register double tmp19_0 = A_01_2 + A_10_2;  
                                     const register double tmp14_1 = A_10_0*w8;  
                                     const register double tmp23_1 = tmp3_0*w14;  
                                     const register double tmp35_1 = A_01_0*w8;  
                                     const register double tmp54_1 = tmp13_0*w8;  
                                     const register double tmp20_1 = tmp9_0*w4;  
                                     const register double tmp25_1 = tmp12_0*w12;  
                                     const register double tmp2_1 = A_01_1*w4;  
                                     const register double tmp44_1 = tmp7_0*w7;  
                                     const register double tmp26_1 = tmp10_0*w4;  
                                     const register double tmp52_1 = tmp18_0*w8;  
                                     const register double tmp48_1 = A_10_1*w7;  
                                     const register double tmp46_1 = A_01_3*w8;  
                                     const register double tmp50_1 = A_01_0*w2;  
                                     const register double tmp8_1 = tmp4_0*w5;  
                                     const register double tmp56_1 = tmp19_0*w8;  
                                     const register double tmp9_1 = tmp2_0*w10;  
                                     const register double tmp19_1 = A_10_3*w2;  
                                     const register double tmp47_1 = A_10_2*w4;  
                                     const register double tmp16_1 = tmp3_0*w0;  
                                     const register double tmp18_1 = tmp1_0*w6;  
                                     const register double tmp31_1 = tmp11_0*w12;  
                                     const register double tmp55_1 = tmp15_0*w2;  
                                     const register double tmp39_1 = A_10_2*w7;  
                                     const register double tmp11_1 = tmp6_0*w7;  
                                     const register double tmp40_1 = tmp11_0*w17;  
                                     const register double tmp34_1 = tmp15_0*w8;  
                                     const register double tmp33_1 = tmp14_0*w5;  
                                     const register double tmp24_1 = tmp11_0*w13;  
                                     const register double tmp3_1 = tmp1_0*w0;  
                                     const register double tmp5_1 = tmp2_0*w3;  
                                     const register double tmp43_1 = tmp17_0*w5;  
                                     const register double tmp15_1 = A_01_2*w4;  
                                     const register double tmp53_1 = tmp19_0*w2;  
                                     const register double tmp27_1 = tmp3_0*w11;  
                                     const register double tmp32_1 = tmp13_0*w2;  
                                     const register double tmp10_1 = tmp5_0*w9;  
                                     const register double tmp37_1 = A_10_1*w4;  
                                     const register double tmp38_1 = tmp5_0*w15;  
                                     const register double tmp17_1 = A_01_1*w7;  
                                     const register double tmp12_1 = tmp7_0*w4;  
                                     const register double tmp22_1 = tmp10_0*w7;  
                                     const register double tmp57_1 = tmp18_0*w2;  
                                     const register double tmp28_1 = tmp9_0*w7;  
                                     const register double tmp29_1 = tmp1_0*w14;  
                                     const register double tmp51_1 = tmp11_0*w16;  
                                     const register double tmp42_1 = tmp12_0*w16;  
                                     const register double tmp49_1 = tmp12_0*w17;  
                                     const register double tmp21_1 = tmp1_0*w11;  
                                     const register double tmp1_1 = tmp0_0*w1;  
                                     const register double tmp45_1 = tmp6_0*w4;  
                                     const register double tmp7_1 = A_10_0*w2;  
                                     const register double tmp6_1 = tmp3_0*w6;  
                                     const register double tmp13_1 = tmp8_0*w1;  
                                     const register double tmp36_1 = tmp16_0*w1;  
                                     const register double tmp41_1 = A_01_3*w2;  
                                     const register double tmp30_1 = tmp12_0*w13;  
                                     const register double tmp4_1 = A_01_2*w7;  
                                     const register double tmp0_1 = A_10_3*w8;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;  
                                 }  
                             }  
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];  
                                     const register double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];  
                                     const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];  
                                     const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];  
                                     const register double tmp0_0 = A_01 + A_10;  
                                     const register double tmp0_1 = A_00*w18;  
                                     const register double tmp10_1 = A_01*w20;  
                                     const register double tmp12_1 = A_00*w26;  
                                     const register double tmp4_1 = A_00*w22;  
                                     const register double tmp8_1 = A_00*w24;  
                                     const register double tmp13_1 = A_10*w19;  
                                     const register double tmp9_1 = tmp0_0*w20;  
                                     const register double tmp3_1 = A_11*w21;  
                                     const register double tmp11_1 = A_11*w27;  
                                     const register double tmp1_1 = A_01*w19;  
                                     const register double tmp6_1 = A_11*w23;  
                                     const register double tmp7_1 = A_11*w25;  
                                     const register double tmp2_1 = A_10*w20;  
                                     const register double tmp5_1 = tmp0_0*w19;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process B //  
                     ///////////////  
                     if (!B.isEmpty()) {  
                         add_EM_S=true;  
                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);  
                         if (B.actsExpanded()) {  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const register double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];  
                                     const register double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];  
                                     const register double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];  
                                     const register double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];  
                                     const register double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];  
                                     const register double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];  
                                     const register double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];  
                                     const register double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];  
                                     const register double tmp3_0 = B_0_0 + B_0_2;  
                                     const register double tmp1_0 = B_1_2 + B_1_3;  
                                     const register double tmp2_0 = B_0_1 + B_0_3;  
                                     const register double tmp0_0 = B_1_0 + B_1_1;  
                                     const register double tmp63_1 = B_1_1*w42;  
                                     const register double tmp79_1 = B_1_1*w40;  
                                     const register double tmp37_1 = tmp3_0*w35;  
                                     const register double tmp8_1 = tmp0_0*w32;  
                                     const register double tmp71_1 = B_0_1*w34;  
                                     const register double tmp19_1 = B_0_3*w31;  
                                     const register double tmp15_1 = B_0_3*w34;  
                                     const register double tmp9_1 = tmp3_0*w34;  
                                     const register double tmp35_1 = B_1_0*w36;  
                                     const register double tmp66_1 = B_0_3*w28;  
                                     const register double tmp28_1 = B_1_0*w42;  
                                     const register double tmp22_1 = B_1_0*w40;  
                                     const register double tmp16_1 = B_1_2*w29;  
                                     const register double tmp6_1 = tmp2_0*w35;  
                                     const register double tmp55_1 = B_1_3*w40;  
                                     const register double tmp50_1 = B_1_3*w42;  
                                     const register double tmp7_1 = tmp1_0*w29;  
                                     const register double tmp1_1 = tmp1_0*w32;  
                                     const register double tmp57_1 = B_0_3*w30;  
                                     const register double tmp18_1 = B_1_1*w32;  
                                     const register double tmp53_1 = B_1_0*w41;  
                                     const register double tmp61_1 = B_1_3*w36;  
                                     const register double tmp27_1 = B_0_3*w38;  
                                     const register double tmp64_1 = B_0_2*w30;  
                                     const register double tmp76_1 = B_0_1*w38;  
                                     const register double tmp39_1 = tmp2_0*w34;  
                                     const register double tmp62_1 = B_0_1*w31;  
                                     const register double tmp56_1 = B_0_0*w31;  
                                     const register double tmp49_1 = B_1_1*w36;  
                                     const register double tmp2_1 = B_0_2*w31;  
                                     const register double tmp23_1 = B_0_2*w33;  
                                     const register double tmp38_1 = B_1_1*w43;  
                                     const register double tmp74_1 = B_1_2*w41;  
                                     const register double tmp43_1 = B_1_1*w41;  
                                     const register double tmp58_1 = B_0_2*w28;  
                                     const register double tmp67_1 = B_0_0*w33;  
                                     const register double tmp33_1 = tmp0_0*w39;  
                                     const register double tmp4_1 = B_0_0*w28;  
                                     const register double tmp20_1 = B_0_0*w30;  
                                     const register double tmp13_1 = B_0_2*w38;  
                                     const register double tmp65_1 = B_1_2*w43;  
                                     const register double tmp0_1 = tmp0_0*w29;  
                                     const register double tmp41_1 = tmp3_0*w33;  
                                     const register double tmp73_1 = B_0_2*w37;  
                                     const register double tmp69_1 = B_0_0*w38;  
                                     const register double tmp48_1 = B_1_2*w39;  
                                     const register double tmp59_1 = B_0_1*w33;  
                                     const register double tmp17_1 = B_1_3*w41;  
                                     const register double tmp5_1 = B_0_3*w33;  
                                     const register double tmp3_1 = B_0_1*w30;  
                                     const register double tmp21_1 = B_0_1*w28;  
                                     const register double tmp42_1 = B_1_0*w29;  
                                     const register double tmp54_1 = B_1_2*w32;  
                                     const register double tmp60_1 = B_1_0*w39;  
                                     const register double tmp32_1 = tmp1_0*w36;  
                                     const register double tmp10_1 = B_0_1*w37;  
                                     const register double tmp14_1 = B_0_0*w35;  
                                     const register double tmp29_1 = B_0_1*w35;  
                                     const register double tmp26_1 = B_1_2*w36;  
                                     const register double tmp30_1 = B_1_3*w43;  
                                     const register double tmp70_1 = B_0_2*w35;  
                                     const register double tmp34_1 = B_1_3*w39;  
                                     const register double tmp51_1 = B_1_0*w43;  
                                     const register double tmp31_1 = B_0_2*w34;  
                                     const register double tmp45_1 = tmp3_0*w28;  
                                     const register double tmp11_1 = tmp1_0*w39;  
                                     const register double tmp52_1 = B_1_1*w29;  
                                     const register double tmp44_1 = B_1_3*w32;  
                                     const register double tmp25_1 = B_1_1*w39;  
                                     const register double tmp47_1 = tmp2_0*w33;  
                                     const register double tmp72_1 = B_1_3*w29;  
                                     const register double tmp40_1 = tmp2_0*w28;  
                                     const register double tmp46_1 = B_1_2*w40;  
                                     const register double tmp36_1 = B_1_2*w42;  
                                     const register double tmp24_1 = B_0_0*w37;  
                                     const register double tmp77_1 = B_0_3*w35;  
                                     const register double tmp68_1 = B_0_3*w37;  
                                     const register double tmp78_1 = B_0_0*w34;  
                                     const register double tmp12_1 = tmp0_0*w36;  
                                     const register double tmp75_1 = B_1_0*w32;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                                 }  
                             }  
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];  
                                     const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];  
                                     const register double tmp6_1 = B_1*w50;  
                                     const register double tmp1_1 = B_1*w45;  
                                     const register double tmp5_1 = B_1*w49;  
                                     const register double tmp4_1 = B_1*w48;  
                                     const register double tmp0_1 = B_0*w44;  
                                     const register double tmp2_1 = B_0*w46;  
                                     const register double tmp7_1 = B_0*w51;  
                                     const register double tmp3_1 = B_0*w47;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;  
                                 }  
                             }  
                         }  
                     }  
                     ///////////////  
                     // process C //  
                     ///////////////  
                     if (!C.isEmpty()) {  
                         add_EM_S=true;  
                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);  
                         if (C.actsExpanded()) {  
                             for (index_t k=0; k<numEq; k++) {  
                                 for (index_t m=0; m<numComp; m++) {  
                                     const register double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];  
                                     const register double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];  
                                     const register double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];  
                                     const register double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];  
                                     const register double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];  
                                     const register double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];  
                                     const register double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];  
                                     const register double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];  
                                     const register double tmp2_0 = C_0_1 + C_0_3;  
                                     const register double tmp1_0 = C_1_2 + C_1_3;  
                                     const register double tmp3_0 = C_0_0 + C_0_2;  
                                     const register double tmp0_0 = C_1_0 + C_1_1;  
                                     const register double tmp64_1 = C_0_2*w30;  
                                     const register double tmp14_1 = C_0_2*w28;  
                                     const register double tmp19_1 = C_0_3*w31;  
                                     const register double tmp22_1 = C_1_0*w40;  
                                     const register double tmp37_1 = tmp3_0*w35;  
                                     const register double tmp29_1 = C_0_1*w35;  
                                     const register double tmp73_1 = C_0_2*w37;  
                                     const register double tmp74_1 = C_1_2*w41;  
                                     const register double tmp52_1 = C_1_3*w39;  
                                     const register double tmp25_1 = C_1_1*w39;  
                                     const register double tmp62_1 = C_0_1*w31;  
                                     const register double tmp79_1 = C_1_1*w40;  
                                     const register double tmp43_1 = C_1_1*w36;  
                                     const register double tmp27_1 = C_0_3*w38;  
                                     const register double tmp28_1 = C_1_0*w42;  
                                     const register double tmp63_1 = C_1_1*w42;  
                                     const register double tmp59_1 = C_0_3*w34;  
                                     const register double tmp72_1 = C_1_3*w29;  
                                     const register double tmp40_1 = tmp2_0*w35;  
                                     const register double tmp13_1 = C_0_3*w30;  
                                     const register double tmp51_1 = C_1_2*w40;  
                                     const register double tmp54_1 = C_1_2*w42;  
                                     const register double tmp12_1 = C_0_0*w31;  
                                     const register double tmp2_1 = tmp1_0*w32;  
                                     const register double tmp68_1 = C_0_2*w31;  
                                     const register double tmp75_1 = C_1_0*w32;  
                                     const register double tmp49_1 = C_1_1*w41;  
                                     const register double tmp4_1 = C_0_2*w35;  
                                     const register double tmp66_1 = C_0_3*w28;  
                                     const register double tmp56_1 = C_0_1*w37;  
                                     const register double tmp5_1 = C_0_1*w34;  
                                     const register double tmp38_1 = tmp2_0*w34;  
                                     const register double tmp76_1 = C_0_1*w38;  
                                     const register double tmp21_1 = C_0_1*w28;  
                                     const register double tmp69_1 = C_0_1*w30;  
                                     const register double tmp53_1 = C_1_0*w36;  
                                     const register double tmp42_1 = C_1_2*w39;  
                                     const register double tmp32_1 = tmp1_0*w29;  
                                     const register double tmp45_1 = C_1_0*w43;  
                                     const register double tmp33_1 = tmp0_0*w32;  
                                     const register double tmp35_1 = C_1_0*w41;  
                                     const register double tmp26_1 = C_1_2*w36;  
                                     const register double tmp67_1 = C_0_0*w33;  
                                     const register double tmp31_1 = C_0_2*w34;  
                                     const register double tmp20_1 = C_0_0*w30;  
                                     const register double tmp70_1 = C_0_0*w28;  
                                     const register double tmp8_1 = tmp0_0*w39;  
                                     const register double tmp30_1 = C_1_3*w43;  
                                     const register double tmp0_1 = tmp0_0*w29;  
                                     const register double tmp17_1 = C_1_3*w41;  
                                     const register double tmp58_1 = C_0_0*w35;  
                                     const register double tmp9_1 = tmp3_0*w33;  
                                     const register double tmp61_1 = C_1_3*w36;  
                                     const register double tmp41_1 = tmp3_0*w34;  
                                     const register double tmp50_1 = C_1_3*w32;  
                                     const register double tmp18_1 = C_1_1*w32;  
                                     const register double tmp6_1 = tmp1_0*w36;  
                                     const register double tmp3_1 = C_0_0*w38;  
                                     const register double tmp34_1 = C_1_1*w29;  
                                     const register double tmp77_1 = C_0_3*w35;  
                                     const register double tmp65_1 = C_1_2*w43;  
                                     const register double tmp71_1 = C_0_3*w33;  
                                     const register double tmp55_1 = C_1_1*w43;  
                                     const register double tmp46_1 = tmp3_0*w28;  
                                     const register double tmp24_1 = C_0_0*w37;  
                                     const register double tmp10_1 = tmp1_0*w39;  
                                     const register double tmp48_1 = C_1_0*w29;  
                                     const register double tmp15_1 = C_0_1*w33;  
                                     const register double tmp36_1 = C_1_2*w32;  
                                     const register double tmp60_1 = C_1_0*w39;  
                                     const register double tmp47_1 = tmp2_0*w33;  
                                     const register double tmp16_1 = C_1_2*w29;  
                                     const register double tmp1_1 = C_0_3*w37;  
                                     const register double tmp7_1 = tmp2_0*w28;  
                                     const register double tmp39_1 = C_1_3*w40;  
                                     const register double tmp44_1 = C_1_3*w42;  
                                     const register double tmp57_1 = C_0_2*w38;  
                                     const register double tmp78_1 = C_0_0*w34;  
                                     const register double tmp11_1 = tmp0_0*w36;  
                                     const register double tmp23_1 = C_0_2*w33;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;