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

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

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

branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp revision 3791 by caltinay, Wed Feb 1 05:10:22 2012 UTC trunk/ripley/src/Rectangle.cpp revision 4375 by caltinay, Mon Apr 22 05:35:52 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2012 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
 extern "C" {  
17  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
18  }  #include <esysUtils/esysFileWriter.h>
19    
20    #ifdef USE_NETCDF
21    #include <netcdfcpp.h>
22    #endif
23    
24  #if USE_SILO  #if USE_SILO
25  #include <silo.h>  #include <silo.h>
# Line 26  extern "C" { Line 31  extern "C" {
31  #include <iomanip>  #include <iomanip>
32    
33  using namespace std;  using namespace std;
34    using esysUtils::FileWriter;
35    
36  namespace ripley {  namespace ripley {
37    
38  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
39                       double y1, int d0, int d1) :                       double y1, int d0, int d1) :
40      RipleyDomain(2),      RipleyDomain(2)
     m_gNE0(n0),  
     m_gNE1(n1),  
     m_x0(x0),  
     m_y0(y0),  
     m_l0(x1-x0),  
     m_l1(y1-y0),  
     m_NX(d0),  
     m_NY(d1)  
41  {  {
42        // ignore subdivision parameters for serial run
43        if (m_mpiInfo->size == 1) {
44            d0=1;
45            d1=1;
46        }
47    
48        bool warn=false;
49        // if number of subdivisions is non-positive, try to subdivide by the same
50        // ratio as the number of elements
51        if (d0<=0 && d1<=0) {
52            warn=true;
53            d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
54            d1=m_mpiInfo->size/d0;
55            if (d0*d1 != m_mpiInfo->size) {
56                // ratios not the same so subdivide side with more elements only
57                if (n0>n1) {
58                    d0=0;
59                    d1=1;
60                } else {
61                    d0=1;
62                    d1=0;
63                }
64            }
65        }
66        if (d0<=0) {
67            // d1 is preset, determine d0 - throw further down if result is no good
68            d0=m_mpiInfo->size/d1;
69        } else if (d1<=0) {
70            // d0 is preset, determine d1 - throw further down if result is no good
71            d1=m_mpiInfo->size/d0;
72        }
73    
74      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
75      // among number of ranks      // among number of ranks
76      if (m_NX*m_NY != m_mpiInfo->size)      if (d0*d1 != m_mpiInfo->size)
77          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
78    
79      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)      if (warn) {
80          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
81                << d1 << "). This may not be optimal!" << endl;
82        }
83    
84        double l0 = x1-x0;
85        double l1 = y1-y0;
86        m_dx[0] = l0/n0;
87        m_dx[1] = l1/n1;
88    
89        if ((n0+1)%d0 > 0) {
90            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
91            l0=m_dx[0]*n0;
92            cout << "Warning: Adjusted number of elements and length. N0="
93                << n0 << ", l0=" << l0 << endl;
94        }
95        if ((n1+1)%d1 > 0) {
96            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
97            l1=m_dx[1]*n1;
98            cout << "Warning: Adjusted number of elements and length. N1="
99                << n1 << ", l1=" << l1 << endl;
100        }
101    
102      if ((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))
103          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
104    
105        m_gNE[0] = n0;
106        m_gNE[1] = n1;
107        m_origin[0] = x0;
108        m_origin[1] = y0;
109        m_length[0] = l0;
110        m_length[1] = l1;
111        m_NX[0] = d0;
112        m_NX[1] = d1;
113    
114      // local number of elements (with and without overlap)      // local number of elements (with and without overlap)
115      m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);      m_NE[0] = m_ownNE[0] = (d0>1 ? (n0+1)/d0 : n0);
116      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)      if (m_mpiInfo->rank%d0>0 && m_mpiInfo->rank%d0<d0-1)
117          m_NE0++;          m_NE[0]++;
118      else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)      else if (d0>1 && m_mpiInfo->rank%d0==d0-1)
119          m_ownNE0--;          m_ownNE[0]--;
120    
121      m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);      m_NE[1] = m_ownNE[1] = (d1>1 ? (n1+1)/d1 : n1);
122      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)      if (m_mpiInfo->rank/d0>0 && m_mpiInfo->rank/d0<d1-1)
123          m_NE1++;          m_NE[1]++;
124      else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)      else if (d1>1 && m_mpiInfo->rank/d0==d1-1)
125          m_ownNE1--;          m_ownNE[1]--;
126    
127      // local number of nodes      // local number of nodes
128      m_N0 = m_NE0+1;      m_NN[0] = m_NE[0]+1;
129      m_N1 = m_NE1+1;      m_NN[1] = m_NE[1]+1;
130    
131      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
132      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);      m_offset[0] = (n0+1)/d0*(m_mpiInfo->rank%d0);
133      if (m_offset0 > 0)      if (m_offset[0] > 0)
134          m_offset0--;          m_offset[0]--;
135      m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);      m_offset[1] = (n1+1)/d1*(m_mpiInfo->rank/d0);
136      if (m_offset1 > 0)      if (m_offset[1] > 0)
137          m_offset1--;          m_offset[1]--;
138    
139      populateSampleIds();      populateSampleIds();
140      createPattern();      createPattern();
# Line 97  bool Rectangle::operator==(const Abstrac Line 156  bool Rectangle::operator==(const Abstrac
156      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
157      if (o) {      if (o) {
158          return (RipleyDomain::operator==(other) &&          return (RipleyDomain::operator==(other) &&
159                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1                  m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1]
160                  && m_x0==o->m_x0 && m_y0==o->m_y0                  && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1]
161                  && m_l0==o->m_l0 && m_l1==o->m_l1                  && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1]
162                  && m_NX==o->m_NX && m_NY==o->m_NY);                  && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1]);
163      }      }
164    
165      return false;      return false;
166  }  }
167    
168    void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
169                const vector<int>& first, const vector<int>& numValues,
170                const vector<int>& multiplier) const
171    {
172    #ifdef USE_NETCDF
173        // check destination function space
174        int myN0, myN1;
175        if (out.getFunctionSpace().getTypeCode() == Nodes) {
176            myN0 = m_NN[0];
177            myN1 = m_NN[1];
178        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
179                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
180            myN0 = m_NE[0];
181            myN1 = m_NE[1];
182        } else
183            throw RipleyException("readNcGrid(): invalid function space for output data object");
184    
185        if (first.size() != 2)
186            throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
187    
188        if (numValues.size() != 2)
189            throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
190    
191        if (multiplier.size() != 2)
192            throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");
193        for (size_t i=0; i<multiplier.size(); i++)
194            if (multiplier[i]<1)
195                throw RipleyException("readNcGrid(): all multipliers must be positive");
196    
197        // check file existence and size
198        NcFile f(filename.c_str(), NcFile::ReadOnly);
199        if (!f.is_valid())
200            throw RipleyException("readNcGrid(): cannot open file");
201    
202        NcVar* var = f.get_var(varname.c_str());
203        if (!var)
204            throw RipleyException("readNcGrid(): invalid variable");
205    
206        // TODO: rank>0 data support
207        const int numComp = out.getDataPointSize();
208        if (numComp > 1)
209            throw RipleyException("readNcGrid(): only scalar data supported");
210    
211        const int dims = var->num_dims();
212        const long *edges = var->edges();
213    
214        // is this a slice of the data object (dims!=2)?
215        // note the expected ordering of edges (as in numpy: y,x)
216        if ( (dims==2 && (numValues[1] > edges[0] || numValues[0] > edges[1]))
217                || (dims==1 && numValues[1]>1) ) {
218            throw RipleyException("readNcGrid(): not enough data in file");
219        }
220    
221        // check if this rank contributes anything
222        if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset[0] ||
223                first[1] >= m_offset[1]+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset[1])
224            return;
225    
226        // now determine how much this rank has to write
227    
228        // first coordinates in data object to write to
229        const int first0 = max(0, first[0]-m_offset[0]);
230        const int first1 = max(0, first[1]-m_offset[1]);
231        // indices to first value in file
232        const int idx0 = max(0, m_offset[0]-first[0]);
233        const int idx1 = max(0, m_offset[1]-first[1]);
234        // number of values to read
235        const int num0 = min(numValues[0]-idx0, myN0-first0);
236        const int num1 = min(numValues[1]-idx1, myN1-first1);
237    
238        vector<double> values(num0*num1);
239        if (dims==2) {
240            var->set_cur(idx1, idx0);
241            var->get(&values[0], num1, num0);
242        } else {
243            var->set_cur(idx0);
244            var->get(&values[0], num0);
245        }
246    
247        const int dpp = out.getNumDataPointsPerSample();
248        out.requireWrite();
249    
250        for (index_t y=0; y<num1; y++) {
251    #pragma omp parallel for
252            for (index_t x=0; x<num0; x++) {
253                const int baseIndex = first0+x*multiplier[0]
254                                      +(first1+y*multiplier[1])*myN0;
255                const int srcIndex = y*num0+x;
256                if (!isnan(values[srcIndex])) {
257                    for (index_t m1=0; m1<multiplier[1]; m1++) {
258                        for (index_t m0=0; m0<multiplier[0]; m0++) {
259                            const int dataIndex = baseIndex+m0+m1*myN0;
260                            double* dest = out.getSampleDataRW(dataIndex);
261                            for (index_t q=0; q<dpp; q++) {
262                                *dest++ = values[srcIndex];
263                            }
264                        }
265                    }
266                }
267            }
268        }
269    #else
270        throw RipleyException("readNcGrid(): not compiled with netCDF support");
271    #endif
272    }
273    
274    void Rectangle::readBinaryGrid(escript::Data& out, string filename,
275                                   const vector<int>& first,
276                                   const vector<int>& numValues,
277                                   const vector<int>& multiplier) const
278    {
279        // check destination function space
280        int myN0, myN1;
281        if (out.getFunctionSpace().getTypeCode() == Nodes) {
282            myN0 = m_NN[0];
283            myN1 = m_NN[1];
284        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
285                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
286            myN0 = m_NE[0];
287            myN1 = m_NE[1];
288        } else
289            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
290    
291        // check file existence and size
292        ifstream f(filename.c_str(), ifstream::binary);
293        if (f.fail()) {
294            throw RipleyException("readBinaryGrid(): cannot open file");
295        }
296        f.seekg(0, ios::end);
297        const int numComp = out.getDataPointSize();
298        const int filesize = f.tellg();
299        const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(float);
300        if (filesize < reqsize) {
301            f.close();
302            throw RipleyException("readBinaryGrid(): not enough data in file");
303        }
304    
305        // check if this rank contributes anything
306        if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0] <= m_offset[0] ||
307                first[1] >= m_offset[1]+myN1 || first[1]+numValues[1] <= m_offset[1]) {
308            f.close();
309            return;
310        }
311    
312        // now determine how much this rank has to write
313    
314        // first coordinates in data object to write to
315        const int first0 = max(0, first[0]-m_offset[0]);
316        const int first1 = max(0, first[1]-m_offset[1]);
317        // indices to first value in file
318        const int idx0 = max(0, m_offset[0]-first[0]);
319        const int idx1 = max(0, m_offset[1]-first[1]);
320        // number of values to read
321        const int num0 = min(numValues[0]-idx0, myN0-first0);
322        const int num1 = min(numValues[1]-idx1, myN1-first1);
323    
324        out.requireWrite();
325        vector<float> values(num0*numComp);
326        const int dpp = out.getNumDataPointsPerSample();
327    
328        for (index_t y=0; y<num1; y++) {
329            const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);
330            f.seekg(fileofs*sizeof(float));
331            f.read((char*)&values[0], num0*numComp*sizeof(float));
332            for (index_t x=0; x<num0; x++) {
333                const int baseIndex = first0+x*multiplier[0]
334                                        +(first1+y*multiplier[1])*myN0;
335                for (index_t m1=0; m1<multiplier[1]; m1++) {
336                    for (index_t m0=0; m0<multiplier[0]; m0++) {
337                        const int dataIndex = baseIndex+m0+m1*myN0;
338                        double* dest = out.getSampleDataRW(dataIndex);
339                        for (index_t c=0; c<numComp; c++) {
340                            if (!std::isnan(values[x*numComp+c])) {
341                                for (index_t q=0; q<dpp; q++) {
342                                    *dest++ = static_cast<double>(values[x*numComp+c]);
343                                }
344                            }
345                        }
346                    }
347                }
348            }
349        }
350    
351        f.close();
352    }
353    
354    void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,
355                                    int byteOrder, int dataType) const
356    {
357        // the mapping is not universally correct but should work on our
358        // supported platforms
359        switch (dataType) {
360            case DATATYPE_INT32:
361                writeBinaryGridImpl<int>(in, filename, byteOrder);
362                break;
363            case DATATYPE_FLOAT32:
364                writeBinaryGridImpl<float>(in, filename, byteOrder);
365                break;
366            case DATATYPE_FLOAT64:
367                writeBinaryGridImpl<double>(in, filename, byteOrder);
368                break;
369            default:
370                throw RipleyException("writeBinaryGrid(): invalid or unsupported datatype");
371        }
372    }
373    
374    template<typename ValueType>
375    void Rectangle::writeBinaryGridImpl(const escript::Data& in,
376                                        const string& filename, int byteOrder) const
377    {
378        // check function space and determine number of points
379        int myN0, myN1;
380        int totalN0, totalN1;
381        if (in.getFunctionSpace().getTypeCode() == Nodes) {
382            myN0 = m_NN[0];
383            myN1 = m_NN[1];
384            totalN0 = m_gNE[0]+1;
385            totalN1 = m_gNE[1]+1;
386        } else if (in.getFunctionSpace().getTypeCode() == Elements ||
387                    in.getFunctionSpace().getTypeCode() == ReducedElements) {
388            myN0 = m_NE[0];
389            myN1 = m_NE[1];
390            totalN0 = m_gNE[0];
391            totalN1 = m_gNE[1];
392        } else
393            throw RipleyException("writeBinaryGrid(): invalid function space of data object");
394    
395        const int numComp = in.getDataPointSize();
396        const int dpp = in.getNumDataPointsPerSample();
397    
398        if (numComp > 1 || dpp > 1)
399            throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
400    
401        escript::Data* _in = const_cast<escript::Data*>(&in);
402        const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;
403    
404        // from here on we know that each sample consists of one value
405        FileWriter* fw = new FileWriter();
406        fw->openFile(filename, fileSize);
407        MPIBarrier();
408    
409        for (index_t y=0; y<myN1; y++) {
410            const int fileofs = (m_offset[0]+(m_offset[1]+y)*totalN0)*sizeof(ValueType);
411            ostringstream oss;
412    
413            for (index_t x=0; x<myN0; x++) {
414                const double* sample = _in->getSampleDataRO(y*myN0+x);
415                ValueType fvalue = static_cast<ValueType>(*sample);
416                if (byteOrder == BYTEORDER_NATIVE) {
417                    oss.write((char*)&fvalue, sizeof(fvalue));
418                } else {
419                    char* value = reinterpret_cast<char*>(&fvalue);
420                    oss.write(byte_swap32(value), sizeof(fvalue));
421                }
422            }
423            fw->writeAt(oss, fileofs);
424        }
425        fw->close();
426    }
427    
428  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
429  {  {
430  #if USE_SILO  #if USE_SILO
# Line 114  void Rectangle::dump(const string& fileN Line 433  void Rectangle::dump(const string& fileN
433          fn+=".silo";          fn+=".silo";
434      }      }
435    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
436      int driver=DB_HDF5;          int driver=DB_HDF5;    
437      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
438        const char* blockDirFmt = "/block%04d";
439    
440  #ifdef ESYS_MPI  #ifdef ESYS_MPI
441      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
442        const int NUM_SILO_FILES = 1;
443  #endif  #endif
444    
445      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 168  void Rectangle::dump(const string& fileN Line 487  void Rectangle::dump(const string& fileN
487      }      }
488      */      */
489    
490      boost::scoped_ptr<double> x(new double[m_N0]);      boost::scoped_ptr<double> x(new double[m_NN[0]]);
491      boost::scoped_ptr<double> y(new double[m_N1]);      boost::scoped_ptr<double> y(new double[m_NN[1]]);
492      double* coords[2] = { x.get(), y.get() };      double* coords[2] = { x.get(), y.get() };
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
493  #pragma omp parallel  #pragma omp parallel
494      {      {
495  #pragma omp for nowait  #pragma omp for nowait
496          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
497              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=getLocalCoordinate(i0, 0);
498          }          }
499  #pragma omp for nowait  #pragma omp for nowait
500          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
501              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=getLocalCoordinate(i1, 1);
502          }          }
503      }      }
504      IndexVector dims = getNumNodesPerDim();      int* dims = const_cast<int*>(getNumNodesPerDim());
505    
506      // write mesh      // write mesh
507      DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,
508              DB_COLLINEAR, NULL);              DB_COLLINEAR, NULL);
509    
510      // write node ids      // write node ids
511      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 2,      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 2,
512              NULL, 0, DB_INT, DB_NODECENT, NULL);              NULL, 0, DB_INT, DB_NODECENT, NULL);
513    
514      // write element ids      // write element ids
515      dims = getNumElementsPerDim();      dims = const_cast<int*>(getNumElementsPerDim());
516      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
517              &dims[0], 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);              dims, 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);
518    
519      // rank 0 writes multimesh and multivar      // rank 0 writes multimesh and multivar
520      if (m_mpiInfo->rank == 0) {      if (m_mpiInfo->rank == 0) {
# Line 289  bool Rectangle::ownSample(int fsType, in Line 606  bool Rectangle::ownSample(int fsType, in
606          case Elements:          case Elements:
607          case ReducedElements:          case ReducedElements:
608              // check ownership of element's bottom left node              // check ownership of element's bottom left node
609              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());
610          case FaceElements:          case FaceElements:
611          case ReducedFaceElements:          case ReducedFaceElements:
612              {              {
613                  // determine which face the sample belongs to before                  // determine which face the sample belongs to before
614                  // checking ownership of corresponding element's first node                  // checking ownership of corresponding element's first node
                 const IndexVector faces = getNumFacesPerBoundary();  
615                  dim_t n=0;                  dim_t n=0;
616                  for (size_t i=0; i<faces.size(); i++) {                  for (size_t i=0; i<4; i++) {
617                      n+=faces[i];                      n+=m_faceCount[i];
618                      if (id<n) {                      if (id<n) {
619                          index_t k;                          index_t k;
620                          if (i==1)                          if (i==1)
621                              k=m_N0-2;                              k=m_NN[0]-2;
622                          else if (i==3)                          else if (i==3)
623                              k=m_N0*(m_N1-2);                              k=m_NN[0]*(m_NN[1]-2);
624                          else                          else
625                              k=0;                              k=0;
626                          // determine whether to move right or up                          // determine whether to move right or up
627                          const index_t delta=(i/2==0 ? m_N0 : 1);                          const index_t delta=(i/2==0 ? m_NN[0] : 1);
628                          return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());                          return (m_dofMap[k+(id-n+m_faceCount[i])*delta] < getNumDOF());
629                      }                      }
630                  }                  }
631                  return false;                  return false;
# Line 331  void Rectangle::setToNormal(escript::Dat Line 647  void Rectangle::setToNormal(escript::Dat
647          {          {
648              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
649  #pragma omp for nowait  #pragma omp for nowait
650                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
651                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
652                      // set vector at two quadrature points                      // set vector at two quadrature points
653                      *o++ = -1.;                      *o++ = -1.;
# Line 343  void Rectangle::setToNormal(escript::Dat Line 659  void Rectangle::setToNormal(escript::Dat
659    
660              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
661  #pragma omp for nowait  #pragma omp for nowait
662                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
663                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
664                      // set vector at two quadrature points                      // set vector at two quadrature points
665                      *o++ = 1.;                      *o++ = 1.;
# Line 355  void Rectangle::setToNormal(escript::Dat Line 671  void Rectangle::setToNormal(escript::Dat
671    
672              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
673  #pragma omp for nowait  #pragma omp for nowait
674                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
675                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
676                      // set vector at two quadrature points                      // set vector at two quadrature points
677                      *o++ = 0.;                      *o++ = 0.;
# Line 367  void Rectangle::setToNormal(escript::Dat Line 683  void Rectangle::setToNormal(escript::Dat
683    
684              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
685  #pragma omp for nowait  #pragma omp for nowait
686                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
687                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
688                      // set vector at two quadrature points                      // set vector at two quadrature points
689                      *o++ = 0.;                      *o++ = 0.;
# Line 383  void Rectangle::setToNormal(escript::Dat Line 699  void Rectangle::setToNormal(escript::Dat
699          {          {
700              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
701  #pragma omp for nowait  #pragma omp for nowait
702                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
703                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
704                      *o++ = -1.;                      *o++ = -1.;
705                      *o = 0.;                      *o = 0.;
# Line 392  void Rectangle::setToNormal(escript::Dat Line 708  void Rectangle::setToNormal(escript::Dat
708    
709              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
710  #pragma omp for nowait  #pragma omp for nowait
711                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
712                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
713                      *o++ = 1.;                      *o++ = 1.;
714                      *o = 0.;                      *o = 0.;
# Line 401  void Rectangle::setToNormal(escript::Dat Line 717  void Rectangle::setToNormal(escript::Dat
717    
718              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
719  #pragma omp for nowait  #pragma omp for nowait
720                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
721                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
722                      *o++ = 0.;                      *o++ = 0.;
723                      *o = -1.;                      *o = -1.;
# Line 410  void Rectangle::setToNormal(escript::Dat Line 726  void Rectangle::setToNormal(escript::Dat
726    
727              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
728  #pragma omp for nowait  #pragma omp for nowait
729                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
730                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
731                      *o++ = 0.;                      *o++ = 0.;
732                      *o = 1.;                      *o = 1.;
# Line 432  void Rectangle::setToSize(escript::Data& Line 748  void Rectangle::setToSize(escript::Data&
748              || out.getFunctionSpace().getTypeCode() == ReducedElements) {              || out.getFunctionSpace().getTypeCode() == ReducedElements) {
749          out.requireWrite();          out.requireWrite();
750          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
751          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);  
752  #pragma omp parallel for  #pragma omp parallel for
753          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
754              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 444  void Rectangle::setToSize(escript::Data& Line 758  void Rectangle::setToSize(escript::Data&
758              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
759          out.requireWrite();          out.requireWrite();
760          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
         const double hSize=getFirstCoordAndSpacing(0).second;  
         const double vSize=getFirstCoordAndSpacing(1).second;  
761  #pragma omp parallel  #pragma omp parallel
762          {          {
763              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
764  #pragma omp for nowait  #pragma omp for nowait
765                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
766                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
767                      fill(o, o+numQuad, vSize);                      fill(o, o+numQuad, m_dx[1]);
768                  }                  }
769              }              }
770    
771              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
772  #pragma omp for nowait  #pragma omp for nowait
773                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
774                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
775                      fill(o, o+numQuad, vSize);                      fill(o, o+numQuad, m_dx[1]);
776                  }                  }
777              }              }
778    
779              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
780  #pragma omp for nowait  #pragma omp for nowait
781                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
782                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
783                      fill(o, o+numQuad, hSize);                      fill(o, o+numQuad, m_dx[0]);
784                  }                  }
785              }              }
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                      fill(o, o+numQuad, hSize);                      fill(o, o+numQuad, m_dx[0]);
792                  }                  }
793              }              }
794          } // end of parallel section          } // end of parallel section
# Line 489  void Rectangle::setToSize(escript::Data& Line 801  void Rectangle::setToSize(escript::Data&
801      }      }
802  }  }
803    
 Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  
                                                 bool reducedColOrder) const  
 {  
     /* FIXME: reduced  
     if (reducedRowOrder || reducedColOrder)  
         throw RipleyException("getPattern() not implemented for reduced order");  
     */  
     return m_pattern;  
 }  
   
804  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
805  {  {
806      RipleyDomain::Print_Mesh_Info(full);      RipleyDomain::Print_Mesh_Info(full);
# Line 506  void Rectangle::Print_Mesh_Info(const bo Line 808  void Rectangle::Print_Mesh_Info(const bo
808          cout << "     Id  Coordinates" << endl;          cout << "     Id  Coordinates" << endl;
809          cout.precision(15);          cout.precision(15);
810          cout.setf(ios::scientific, ios::floatfield);          cout.setf(ios::scientific, ios::floatfield);
         pair<double,double> xdx = getFirstCoordAndSpacing(0);  
         pair<double,double> ydy = getFirstCoordAndSpacing(1);  
811          for (index_t i=0; i < getNumNodes(); i++) {          for (index_t i=0; i < getNumNodes(); i++) {
812              cout << "  " << setw(5) << m_nodeId[i]              cout << "  " << setw(5) << m_nodeId[i]
813                  << "  " << xdx.first+(i%m_N0)*xdx.second                  << "  " << getLocalCoordinate(i%m_NN[0], 0)
814                  << "  " << ydy.first+(i/m_N0)*ydy.second << endl;                  << "  " << getLocalCoordinate(i/m_NN[0], 1) << endl;
815          }          }
816      }      }
817  }  }
818    
 IndexVector Rectangle::getNumNodesPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_N0);  
     ret.push_back(m_N1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumElementsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NE0);  
     ret.push_back(m_NE1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumFacesPerBoundary() const  
 {  
     IndexVector ret(4, 0);  
     //left  
     if (m_offset0==0)  
         ret[0]=m_NE1;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         ret[1]=m_NE1;  
     //bottom  
     if (m_offset1==0)  
         ret[2]=m_NE0;  
     //top  
     if (m_mpiInfo->rank/m_NX==m_NY-1)  
         ret[3]=m_NE0;  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumSubdivisionsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NX);  
     ret.push_back(m_NY);  
     return ret;  
 }  
   
 pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  
 {  
     if (dim==0) {  
         return pair<double,double>(m_x0+(m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);  
     } else if (dim==1) {  
         return pair<double,double>(m_y0+(m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);  
     }  
     throw RipleyException("getFirstCoordAndSpacing: invalid argument");  
 }  
   
 //protected  
 dim_t Rectangle::getNumDOF() const  
 {  
     return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;  
 }  
   
 //protected  
 dim_t Rectangle::getNumFaceElements() const  
 {  
     const IndexVector faces = getNumFacesPerBoundary();  
     dim_t n=0;  
     for (size_t i=0; i<faces.size(); i++)  
         n+=faces[i];  
     return n;  
 }  
819    
820  //protected  //protected
821  void Rectangle::assembleCoordinates(escript::Data& arg) const  void Rectangle::assembleCoordinates(escript::Data& arg) const
# Line 594  void Rectangle::assembleCoordinates(escr Line 827  void Rectangle::assembleCoordinates(escr
827      if (!numSamplesEqual(&x, 1, getNumNodes()))      if (!numSamplesEqual(&x, 1, getNumNodes()))
828          throw RipleyException("setToX: Illegal number of samples in Data object");          throw RipleyException("setToX: Illegal number of samples in Data object");
829    
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
830      arg.requireWrite();      arg.requireWrite();
831  #pragma omp parallel for  #pragma omp parallel for
832      for (dim_t i1 = 0; i1 < m_N1; i1++) {      for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
833          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
834              double* point = arg.getSampleDataRW(i0+m_N0*i1);              double* point = arg.getSampleDataRW(i0+m_NN[0]*i1);
835              point[0] = xdx.first+i0*xdx.second;              point[0] = getLocalCoordinate(i0, 0);
836              point[1] = ydy.first+i1*ydy.second;              point[1] = getLocalCoordinate(i1, 1);
837          }          }
838      }      }
839  }  }
# Line 611  void Rectangle::assembleCoordinates(escr Line 842  void Rectangle::assembleCoordinates(escr
842  void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const  void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
843  {  {
844      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
845      const double h0 = m_l0/m_gNE0;      const double cx0 = .21132486540518711775/m_dx[0];
846      const double h1 = m_l1/m_gNE1;      const double cx1 = .78867513459481288225/m_dx[0];
847      const double cx0 = -1./h0;      const double cx2 = 1./m_dx[0];
848      const double cx1 = -.78867513459481288225/h0;      const double cy0 = .21132486540518711775/m_dx[1];
849      const double cx2 = -.5/h0;      const double cy1 = .78867513459481288225/m_dx[1];
850      const double cx3 = -.21132486540518711775/h0;      const double cy2 = 1./m_dx[1];
     const double cx4 = .21132486540518711775/h0;  
     const double cx5 = .5/h0;  
     const double cx6 = .78867513459481288225/h0;  
     const double cx7 = 1./h0;  
     const double cy0 = -1./h1;  
     const double cy1 = -.78867513459481288225/h1;  
     const double cy2 = -.5/h1;  
     const double cy3 = -.21132486540518711775/h1;  
     const double cy4 = .21132486540518711775/h1;  
     const double cy5 = .5/h1;  
     const double cy6 = .78867513459481288225/h1;  
     const double cy7 = 1./h1;  
851    
852      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
853          out.requireWrite();          out.requireWrite();
854  #pragma omp parallel for  #pragma omp parallel
855          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
856              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
857                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
858                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
859                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
860                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
861                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
862                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
863                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
864                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
865                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
866                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
867                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
868                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      for (index_t i=0; i < numComp; ++i) {
869                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
870                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
871                  } /* end of component loop i */                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
872              } /* end of k0 loop */                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
873          } /* end of k1 loop */                          o[INDEX3(i,0,2,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
874                            o[INDEX3(i,1,2,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
875                            o[INDEX3(i,0,3,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
876                            o[INDEX3(i,1,3,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
877                        } // end of component loop i
878                    } // end of k0 loop
879                } // end of k1 loop
880            } // end of parallel section
881      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
882          out.requireWrite();          out.requireWrite();
883  #pragma omp parallel for  #pragma omp parallel
884          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
885              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
886                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
887                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
888                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
889                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
890                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
891                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
892                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
893                      o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
894                  } /* end of component loop i */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
895              } /* end of k0 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
896          } /* end of k1 loop */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
897                        for (index_t i=0; i < numComp; ++i) {
898                            o[INDEX3(i,0,0,numComp,2)] = (f_10[i] + f_11[i] - f_00[i] - f_01[i])*cx2/2;
899                            o[INDEX3(i,1,0,numComp,2)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
900                        } // end of component loop i
901                    } // end of k0 loop
902                } // end of k1 loop
903            } // end of parallel section
904      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
905          out.requireWrite();          out.requireWrite();
906  #pragma omp parallel  #pragma omp parallel
907          {          {
908                vector<double> f_00(numComp);
909                vector<double> f_01(numComp);
910                vector<double> f_10(numComp);
911                vector<double> f_11(numComp);
912              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
913  #pragma omp for nowait  #pragma omp for nowait
914                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
915                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
916                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
917                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
918                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
919                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
920                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
921                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
922                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy2;
923                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
924                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy2;
925                      } /* end of component loop i */                      } // end of component loop i
926                  } /* end of k1 loop */                  } // end of k1 loop
927              } /* end of face 0 */              } // end of face 0
928              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
929  #pragma omp for nowait  #pragma omp for nowait
930                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
931                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
932                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
933                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
934                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
935                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
936                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
937                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
938                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_11[i]-f_10[i])*cy2;
939                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
940                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = (f_11[i]-f_10[i])*cy2;
941                      } /* end of component loop i */                      } // end of component loop i
942                  } /* end of k1 loop */                  } // end of k1 loop
943              } /* end of face 1 */              } // end of face 1
944              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
945  #pragma omp for nowait  #pragma omp for nowait
946                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
947                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
948                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
949                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
950                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
951                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
952                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
953                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx2;
954                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
955                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx2;
956                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
957                      } /* end of component loop i */                      } // end of component loop i
958                  } /* end of k0 loop */                  } // end of k0 loop
959              } /* end of face 2 */              } // end of face 2
960              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
961  #pragma omp for nowait  #pragma omp for nowait
962                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
963                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
964                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
965                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
966                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
967                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
968                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
969                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_11[i]-f_01[i])*cx2;
970                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
971                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = (f_11[i]-f_01[i])*cx2;
972                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
973                      } /* end of component loop i */                      } // end of component loop i
974                  } /* end of k0 loop */                  } // end of k0 loop
975              } /* end of face 3 */              } // end of face 3
976          } // end of parallel section          } // end of parallel section
977    
978      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
979          out.requireWrite();          out.requireWrite();
980  #pragma omp parallel  #pragma omp parallel
981          {          {
982                vector<double> f_00(numComp);
983                vector<double> f_01(numComp);
984                vector<double> f_10(numComp);
985                vector<double> f_11(numComp);
986              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
987  #pragma omp for nowait  #pragma omp for nowait
988                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
989                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
990                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
991                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
992                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
993                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
994                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
995                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i] + f_11[i] - f_00[i] - f_01[i])*cx2/2;
996                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy2;
997                      } /* end of component loop i */                      } // end of component loop i
998                  } /* end of k1 loop */                  } // end of k1 loop
999              } /* end of face 0 */              } // end of face 0
1000              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1001  #pragma omp for nowait  #pragma omp for nowait
1002                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1003                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
1004                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
1005                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1006                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1007                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1008                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1009                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i] + f_11[i] - f_00[i] - f_01[i])*cx2/2;
1010                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_11[i]-f_10[i])*cy2;
1011                      } /* end of component loop i */                      } // end of component loop i
1012                  } /* end of k1 loop */                  } // end of k1 loop
1013              } /* end of face 1 */              } // end of face 1
1014              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1015  #pragma omp for nowait  #pragma omp for nowait
1016                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1017                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1018                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
1019                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1020                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
1021                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1022                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1023                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx2;
1024                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
1025                      } /* end of component loop i */                      } // end of component loop i
1026                  } /* end of k0 loop */                  } // end of k0 loop
1027              } /* end of face 2 */              } // end of face 2
1028              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1029  #pragma omp for nowait  #pragma omp for nowait
1030                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1031                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1032                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1033                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1034                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1035                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1036                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1037                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_11[i]-f_01[i])*cx2;
1038                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
1039                      } /* end of component loop i */                      } // end of component loop i
1040                  } /* end of k0 loop */                  } // end of k0 loop
1041              } /* end of face 3 */              } // end of face 3
1042          } // end of parallel section          } // end of parallel section
1043      }      }
1044  }  }
# Line 805  void Rectangle::assembleGradient(escript Line 1047  void Rectangle::assembleGradient(escript
1047  void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const  void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1048  {  {
1049      const dim_t numComp = arg.getDataPointSize();      const dim_t numComp = arg.getDataPointSize();
1050      const double h0 = m_l0/m_gNE0;      const index_t left = (m_offset[0]==0 ? 0 : 1);
1051      const double h1 = m_l1/m_gNE1;      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1052      const index_t left = (m_offset0==0 ? 0 : 1);      const int fs=arg.getFunctionSpace().getTypeCode();
1053      const index_t bottom = (m_offset1==0 ? 0 : 1);      if (fs == Elements && arg.actsExpanded()) {
     if (arg.getFunctionSpace().getTypeCode() == Elements) {  
         const double w = h0*h1/4.;  
1054  #pragma omp parallel  #pragma omp parallel
1055          {          {
1056              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1057                const double w = m_dx[0]*m_dx[1]/4.;
1058  #pragma omp for nowait  #pragma omp for nowait
1059              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1060                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1061                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1062                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1063                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1064                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1065                          const double f2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
1066                          const double f3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
1067                          int_local[i]+=(f0+f1+f2+f3)*w;                          int_local[i]+=(f0+f1+f2+f3)*w;
1068                      }  /* end of component loop i */                      }  // end of component loop i
1069                  } /* end of k0 loop */                  } // end of k0 loop
1070              } /* end of k1 loop */              } // end of k1 loop
   
1071  #pragma omp critical  #pragma omp critical
1072              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1073                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1074          } // end of parallel section          } // end of parallel section
1075      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
1076          const double w = h0*h1;      } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1077            const double w = m_dx[0]*m_dx[1];
1078  #pragma omp parallel  #pragma omp parallel
1079          {          {
1080              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1081  #pragma omp for nowait  #pragma omp for nowait
1082              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1083                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1084                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1085                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1086                          int_local[i]+=f[i]*w;                          int_local[i]+=f[i]*w;
1087                      }  /* end of component loop i */                      }
1088                  } /* end of k0 loop */                  }
1089              } /* end of k1 loop */              }
   
1090  #pragma omp critical  #pragma omp critical
1091              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1092                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1093          } // end of parallel section          } // end of parallel section
1094      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
1095          const double w0 = h0/2.;      } else if (fs == FaceElements && arg.actsExpanded()) {
         const double w1 = h1/2.;  
1096  #pragma omp parallel  #pragma omp parallel
1097          {          {
1098              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1099                const double w0 = m_dx[0]/2.;
1100                const double w1 = m_dx[1]/2.;
1101              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1102  #pragma omp for nowait  #pragma omp for nowait
1103                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1104                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1105                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1106                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1107                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1108                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1109                      }  /* end of component loop i */                      }  // end of component loop i
1110                  } /* end of k1 loop */                  } // end of k1 loop
1111              }              }
1112    
1113              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1114  #pragma omp for nowait  #pragma omp for nowait
1115                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1116                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1117                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1118                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1119                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1120                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1121                      }  /* end of component loop i */                      }  // end of component loop i
1122                  } /* end of k1 loop */                  } // end of k1 loop
1123              }              }
1124    
1125              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1126  #pragma omp for nowait  #pragma omp for nowait
1127                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1128                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1129                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1130                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1131                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1132                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1133                      }  /* end of component loop i */                      }  // end of component loop i
1134                  } /* end of k0 loop */                  } // end of k0 loop
1135              }              }
1136    
1137              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1138  #pragma omp for nowait  #pragma omp for nowait
1139                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1140                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1141                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1142                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1143                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1144                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1145                      }  /* end of component loop i */                      }  // end of component loop i
1146                  } /* end of k0 loop */                  } // end of k0 loop
1147              }              }
   
1148  #pragma omp critical  #pragma omp critical
1149              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1150                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1151          } // end of parallel section          } // end of parallel section
1152      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
1153        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1154  #pragma omp parallel  #pragma omp parallel
1155          {          {
1156              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1157              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1158  #pragma omp for nowait  #pragma omp for nowait
1159                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1160                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1161                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1162                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*m_dx[1];
1163                      }  /* end of component loop i */                      }
1164                  } /* end of k1 loop */                  }
1165              }              }
1166    
1167              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1168  #pragma omp for nowait  #pragma omp for nowait
1169                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1170                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1171                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1172                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*m_dx[1];
1173                      }  /* end of component loop i */                      }
1174                  } /* end of k1 loop */                  }
1175              }              }
1176    
1177              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1178  #pragma omp for nowait  #pragma omp for nowait
1179                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1180                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1181                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1182                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*m_dx[0];
1183                      }  /* end of component loop i */                      }
1184                  } /* end of k0 loop */                  }
1185              }              }
1186    
1187              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1188  #pragma omp for nowait  #pragma omp for nowait
1189                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1190                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1191                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1192                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*m_dx[0];
1193                      }  /* end of component loop i */                      }
1194                  } /* end of k0 loop */                  }
1195              }              }
1196    
1197  #pragma omp critical  #pragma omp critical
1198              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1199                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1200          } // end of parallel section          } // end of parallel section
1201      }      } // function space selector
1202  }  }
1203    
1204  //protected  //protected
1205  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
1206  {  {
1207      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1208      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1209      const int x=node%nDOF0;      const int x=node%nDOF0;
1210      const int y=node/nDOF0;      const int y=node/nDOF0;
1211      dim_t num=0;      dim_t num=0;
# Line 994  void Rectangle::nodesToDOF(escript::Data Line 1235  void Rectangle::nodesToDOF(escript::Data
1235      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1236      out.requireWrite();      out.requireWrite();
1237    
1238      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1239      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1240      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1241      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1242  #pragma omp parallel for  #pragma omp parallel for
1243      for (index_t i=0; i<nDOF1; i++) {      for (index_t i=0; i<nDOF1; i++) {
1244          for (index_t j=0; j<nDOF0; j++) {          for (index_t j=0; j<nDOF0; j++) {
1245              const index_t n=j+left+(i+bottom)*m_N0;              const index_t n=j+left+(i+bottom)*m_NN[0];
1246              const double* src=in.getSampleDataRO(n);              const double* src=in.getSampleDataRO(n);
1247              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
1248          }          }
# Line 1027  void Rectangle::dofToNodes(escript::Data Line 1268  void Rectangle::dofToNodes(escript::Data
1268                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1269          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1270      }      }
1271        Paso_Coupler_free(coupler);
1272  }  }
1273    
1274  //private  //private
1275  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
1276  {  {
1277      // identifiers are ordered from left to right, bottom to top globablly.      // degrees of freedom are numbered from left to right, bottom to top in
1278        // each rank, continuing on the next rank (ranks also go left-right,
1279        // bottom-top).
1280        // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
1281        // helps when writing out data rank after rank.
1282    
1283      // build node distribution vector first.      // build node distribution vector first.
1284      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
1285        // constant for all ranks in this implementation
1286      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1287      const dim_t numDOF=getNumDOF();      const dim_t numDOF=getNumDOF();
1288      for (dim_t k=1; k<m_mpiInfo->size; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
# Line 1045  void Rectangle::populateSampleIds() Line 1292  void Rectangle::populateSampleIds()
1292      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
1293      m_dofId.resize(numDOF);      m_dofId.resize(numDOF);
1294      m_elementId.resize(getNumElements());      m_elementId.resize(getNumElements());
1295    
1296        // populate face element counts
1297        //left
1298        if (m_offset[0]==0)
1299            m_faceCount[0]=m_NE[1];
1300        else
1301            m_faceCount[0]=0;
1302        //right
1303        if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
1304            m_faceCount[1]=m_NE[1];
1305        else
1306            m_faceCount[1]=0;
1307        //bottom
1308        if (m_offset[1]==0)
1309            m_faceCount[2]=m_NE[0];
1310        else
1311            m_faceCount[2]=0;
1312        //top
1313        if (m_mpiInfo->rank/m_NX[0]==m_NX[1]-1)
1314            m_faceCount[3]=m_NE[0];
1315        else
1316            m_faceCount[3]=0;
1317    
1318      m_faceId.resize(getNumFaceElements());      m_faceId.resize(getNumFaceElements());
1319    
1320        const index_t left = (m_offset[0]==0 ? 0 : 1);
1321        const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1322        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1323        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1324    
1325    #define globalNodeId(x,y) \
1326        ((m_offset[0]+x)/nDOF0)*nDOF0*nDOF1+(m_offset[0]+x)%nDOF0 \
1327        + ((m_offset[1]+y)/nDOF1)*nDOF0*nDOF1*m_NX[0]+((m_offset[1]+y)%nDOF1)*nDOF0
1328    
1329        // set corner id's outside the parallel region
1330        m_nodeId[0] = globalNodeId(0, 0);
1331        m_nodeId[m_NN[0]-1] = globalNodeId(m_NN[0]-1, 0);
1332        m_nodeId[m_NN[0]*(m_NN[1]-1)] = globalNodeId(0, m_NN[1]-1);
1333        m_nodeId[m_NN[0]*m_NN[1]-1] = globalNodeId(m_NN[0]-1,m_NN[1]-1);
1334    #undef globalNodeId
1335    
1336  #pragma omp parallel  #pragma omp parallel
1337      {      {
1338          // nodes          // populate degrees of freedom and own nodes (identical id)
1339  #pragma omp for nowait  #pragma omp for nowait
1340          for (dim_t i1=0; i1<m_N1; i1++) {          for (dim_t i=0; i<nDOF1; i++) {
1341              for (dim_t i0=0; i0<m_N0; i0++) {              for (dim_t j=0; j<nDOF0; j++) {
1342                  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];
1343                    const index_t dofIdx=j+i*nDOF0;
1344                    m_dofId[dofIdx] = m_nodeId[nodeIdx]
1345                        = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
1346              }              }
1347          }          }
1348    
1349          // degrees of freedom          // populate the rest of the nodes (shared with other ranks)
1350            if (m_faceCount[0]==0) { // left column
1351    #pragma omp for nowait
1352                for (dim_t i=0; i<nDOF1; i++) {
1353                    const index_t nodeIdx=(i+bottom)*m_NN[0];
1354                    const index_t dofId=(i+1)*nDOF0-1;
1355                    m_nodeId[nodeIdx]
1356                        = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
1357                }
1358            }
1359            if (m_faceCount[1]==0) { // right column
1360  #pragma omp for nowait  #pragma omp for nowait
1361          for (dim_t k=0; k<numDOF; k++)              for (dim_t i=0; i<nDOF1; i++) {
1362              m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;                  const index_t nodeIdx=(i+bottom+1)*m_NN[0]-1;
1363                    const index_t dofId=i*nDOF0;
1364                    m_nodeId[nodeIdx]
1365                        = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
1366                }
1367            }
1368            if (m_faceCount[2]==0) { // bottom row
1369    #pragma omp for nowait
1370                for (dim_t i=0; i<nDOF0; i++) {
1371                    const index_t nodeIdx=i+left;
1372                    const index_t dofId=nDOF0*(nDOF1-1)+i;
1373                    m_nodeId[nodeIdx]
1374                        = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
1375                }
1376            }
1377            if (m_faceCount[3]==0) { // top row
1378    #pragma omp for nowait
1379                for (dim_t i=0; i<nDOF0; i++) {
1380                    const index_t nodeIdx=m_NN[0]*(m_NN[1]-1)+i+left;
1381                    const index_t dofId=i;
1382                    m_nodeId[nodeIdx]
1383                        = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
1384                }
1385            }
1386    
1387          // elements          // populate element id's
1388  #pragma omp for nowait  #pragma omp for nowait
1389          for (dim_t i1=0; i1<m_NE1; i1++) {          for (dim_t i1=0; i1<m_NE[1]; i1++) {
1390              for (dim_t i0=0; i0<m_NE0; i0++) {              for (dim_t i0=0; i0<m_NE[0]; i0++) {
1391                  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;
1392              }              }
1393          }          }
1394    
# Line 1083  void Rectangle::populateSampleIds() Line 1405  void Rectangle::populateSampleIds()
1405      updateTagsInUse(Elements);      updateTagsInUse(Elements);
1406    
1407      // generate face offset vector and set face tags      // generate face offset vector and set face tags
     const IndexVector facesPerEdge = getNumFacesPerBoundary();  
1408      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
1409      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
1410      m_faceOffset.assign(facesPerEdge.size(), -1);      m_faceOffset.assign(4, -1);
1411      m_faceTags.clear();      m_faceTags.clear();
1412      index_t offset=0;      index_t offset=0;
1413      for (size_t i=0; i<facesPerEdge.size(); i++) {      for (size_t i=0; i<4; i++) {
1414          if (facesPerEdge[i]>0) {          if (m_faceCount[i]>0) {
1415              m_faceOffset[i]=offset;              m_faceOffset[i]=offset;
1416              offset+=facesPerEdge[i];              offset+=m_faceCount[i];
1417              m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);              m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
1418          }          }
1419      }      }
1420      setTagMap("left", LEFT);      setTagMap("left", LEFT);
# Line 1106  void Rectangle::populateSampleIds() Line 1427  void Rectangle::populateSampleIds()
1427  //private  //private
1428  void Rectangle::createPattern()  void Rectangle::createPattern()
1429  {  {
1430      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1431      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1432      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1433      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1434    
1435      // populate node->DOF mapping with own degrees of freedom.      // populate node->DOF mapping with own degrees of freedom.
1436      // The rest is assigned in the loop further down      // The rest is assigned in the loop further down
# Line 1117  void Rectangle::createPattern() Line 1438  void Rectangle::createPattern()
1438  #pragma omp parallel for  #pragma omp parallel for
1439      for (index_t i=bottom; i<bottom+nDOF1; i++) {      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1440          for (index_t j=left; j<left+nDOF0; j++) {          for (index_t j=left; j<left+nDOF0; j++) {
1441              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;              m_dofMap[i*m_NN[0]+j]=(i-bottom)*nDOF0+j-left;
1442          }          }
1443      }      }
1444    
# Line 1130  void Rectangle::createPattern() Line 1451  void Rectangle::createPattern()
1451      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
1452      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
1453      int numShared=0;      int numShared=0;
1454      const int x=m_mpiInfo->rank%m_NX;      const int x=m_mpiInfo->rank%m_NX[0];
1455      const int y=m_mpiInfo->rank/m_NX;      const int y=m_mpiInfo->rank/m_NX[0];
1456      for (int i1=-1; i1<2; i1++) {      for (int i1=-1; i1<2; i1++) {
1457          for (int i0=-1; i0<2; i0++) {          for (int i0=-1; i0<2; i0++) {
1458              // skip this rank              // skip this rank
# Line 1140  void Rectangle::createPattern() Line 1461  void Rectangle::createPattern()
1461              // location of neighbour rank              // location of neighbour rank
1462              const int nx=x+i0;              const int nx=x+i0;
1463              const int ny=y+i1;              const int ny=y+i1;
1464              if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {              if (nx>=0 && ny>=0 && nx<m_NX[0] && ny<m_NX[1]) {
1465                  neighbour.push_back(ny*m_NX+nx);                  neighbour.push_back(ny*m_NX[0]+nx);
1466                  if (i0==0) {                  if (i0==0) {
1467                      // sharing top or bottom edge                      // sharing top or bottom edge
1468                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1469                      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);
1470                      offsetInShared.push_back(offsetInShared.back()+nDOF0);                      offsetInShared.push_back(offsetInShared.back()+nDOF0);
1471                      for (dim_t i=0; i<nDOF0; i++, numShared++) {                      for (dim_t i=0; i<nDOF0; i++, numShared++) {
1472                          sendShared.push_back(firstDOF+i);                          sendShared.push_back(firstDOF+i);
# Line 1160  void Rectangle::createPattern() Line 1481  void Rectangle::createPattern()
1481                  } else if (i1==0) {                  } else if (i1==0) {
1482                      // sharing left or right edge                      // sharing left or right edge
1483                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1484                      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);
1485                      offsetInShared.push_back(offsetInShared.back()+nDOF1);                      offsetInShared.push_back(offsetInShared.back()+nDOF1);
1486                      for (dim_t i=0; i<nDOF1; i++, numShared++) {                      for (dim_t i=0; i<nDOF1; i++, numShared++) {
1487                          sendShared.push_back(firstDOF+i*nDOF0);                          sendShared.push_back(firstDOF+i*nDOF0);
# Line 1170  void Rectangle::createPattern() Line 1491  void Rectangle::createPattern()
1491                          colIndices[firstDOF+i*nDOF0].push_back(numShared);                          colIndices[firstDOF+i*nDOF0].push_back(numShared);
1492                          if (i<nDOF1-1)                          if (i<nDOF1-1)
1493                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1494                          m_dofMap[firstNode+i*m_N0]=numDOF+numShared;                          m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
1495                      }                      }
1496                  } else {                  } else {
1497                      // sharing a node                      // sharing a node
1498                      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);
1499                      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);
1500                      offsetInShared.push_back(offsetInShared.back()+1);                      offsetInShared.push_back(offsetInShared.back()+1);
1501                      sendShared.push_back(dof);                      sendShared.push_back(dof);
1502                      recvShared.push_back(numDOF+numShared);                      recvShared.push_back(numDOF+numShared);
# Line 1285  void Rectangle::addToMatrixAndRHS(Paso_S Line 1606  void Rectangle::addToMatrixAndRHS(Paso_S
1606      IndexVector rowIndex;      IndexVector rowIndex;
1607      rowIndex.push_back(m_dofMap[firstNode]);      rowIndex.push_back(m_dofMap[firstNode]);
1608      rowIndex.push_back(m_dofMap[firstNode+1]);      rowIndex.push_back(m_dofMap[firstNode+1]);
1609      rowIndex.push_back(m_dofMap[firstNode+m_N0]);      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]]);
1610      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]+1]);
1611      if (addF) {      if (addF) {
1612          double *F_p=F.getSampleDataRW(0);          double *F_p=F.getSampleDataRW(0);
1613          for (index_t i=0; i<rowIndex.size(); i++) {          for (index_t i=0; i<rowIndex.size(); i++) {
# Line 1309  void Rectangle::interpolateNodesOnElemen Line 1630  void Rectangle::interpolateNodesOnElemen
1630      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1631      if (reduced) {      if (reduced) {
1632          out.requireWrite();          out.requireWrite();
1633          const double c0 = .25;          const double c0 = 0.25;
1634  #pragma omp parallel for  #pragma omp parallel
1635          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1636              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1637                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_01(numComp);
1638                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1639                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              vector<double> f_11(numComp);
1640                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1641                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1642                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1643                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
1644                  } /* end of component loop i */                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
1645              } /* end of k0 loop */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
1646          } /* end of k1 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
1647                        double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1648                        for (index_t i=0; i < numComp; ++i) {
1649                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1650                        } /* end of component loop i */
1651                    } /* end of k0 loop */
1652                } /* end of k1 loop */
1653            } /* end of parallel section */
1654      } else {      } else {
1655          out.requireWrite();          out.requireWrite();
1656          const double c0 = .16666666666666666667;          const double c0 = 0.16666666666666666667;
1657          const double c1 = .044658198738520451079;          const double c1 = 0.044658198738520451079;
1658          const double c2 = .62200846792814621559;          const double c2 = 0.62200846792814621559;
1659  #pragma omp parallel for  #pragma omp parallel
1660          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1661              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1662                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
1663                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1664                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
1665                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1666                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1667                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1668                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
1669                      o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
1670                      o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
1671                      o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
1672                  } /* end of component loop i */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1673              } /* end of k0 loop */                      for (index_t i=0; i < numComp; ++i) {
1674          } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1675                            o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1676                            o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1677                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1678                        } /* end of component loop i */
1679                    } /* end of k0 loop */
1680                } /* end of k1 loop */
1681            } /* end of parallel section */
1682      }      }
1683  }  }
1684    
# Line 1354  void Rectangle::interpolateNodesOnFaces( Line 1689  void Rectangle::interpolateNodesOnFaces(
1689      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1690      if (reduced) {      if (reduced) {
1691          out.requireWrite();          out.requireWrite();
         const double c0 = .5;  
1692  #pragma omp parallel  #pragma omp parallel
1693          {          {
1694                vector<double> f_00(numComp);
1695                vector<double> f_01(numComp);
1696                vector<double> f_10(numComp);
1697                vector<double> f_11(numComp);
1698              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1699  #pragma omp for nowait  #pragma omp for nowait
1700                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1701                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1702                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1703                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1704                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1705                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = (f_00[i] + f_01[i])/2;
1706                      } /* end of component loop i */                      } /* end of component loop i */
1707                  } /* end of k1 loop */                  } /* end of k1 loop */
1708              } /* end of face 0 */              } /* end of face 0 */
1709              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1710  #pragma omp for nowait  #pragma omp for nowait
1711                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1712                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1713                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1714                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1715                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1716                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = (f_10[i] + f_11[i])/2;
1717                      } /* end of component loop i */                      } /* end of component loop i */
1718                  } /* end of k1 loop */                  } /* end of k1 loop */
1719              } /* end of face 1 */              } /* end of face 1 */
1720              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1721  #pragma omp for nowait  #pragma omp for nowait
1722                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1723                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1724                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1725                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1726                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1727                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = (f_00[i] + f_10[i])/2;
1728                      } /* end of component loop i */                      } /* end of component loop i */
1729                  } /* end of k0 loop */                  } /* end of k0 loop */
1730              } /* end of face 2 */              } /* end of face 2 */
1731              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1732  #pragma omp for nowait  #pragma omp for nowait
1733                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1734                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1735                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1736                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1737                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1738                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = (f_01[i] + f_11[i])/2;
1739                      } /* end of component loop i */                      } /* end of component loop i */
1740                  } /* end of k0 loop */                  } /* end of k0 loop */
1741              } /* end of face 3 */              } /* end of face 3 */
1742          } // end of parallel section          } /* end of parallel section */
1743      } else {      } else {
1744          out.requireWrite();          out.requireWrite();
1745          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1746          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1747  #pragma omp parallel  #pragma omp parallel
1748          {          {
1749                vector<double> f_00(numComp);
1750                vector<double> f_01(numComp);
1751                vector<double> f_10(numComp);
1752                vector<double> f_11(numComp);
1753              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1754  #pragma omp for nowait  #pragma omp for nowait
1755                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1756                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1757                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1758                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1759                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1760                          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];
1761                          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];
1762                      } /* end of component loop i */                      } /* end of component loop i */
1763                  } /* end of k1 loop */                  } /* end of k1 loop */
1764              } /* end of face 0 */              } /* end of face 0 */
1765              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1766  #pragma omp for nowait  #pragma omp for nowait
1767                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1768                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1769                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1770                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1771                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1772                          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];
1773                          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];
1774                      } /* end of component loop i */                      } /* end of component loop i */
1775                  } /* end of k1 loop */                  } /* end of k1 loop */
1776              } /* end of face 1 */              } /* end of face 1 */
1777              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1778  #pragma omp for nowait  #pragma omp for nowait
1779                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1780                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1781                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1782                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1783                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1784                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_10[i] + c1*f_00[i];
1785                          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];
1786                      } /* end of component loop i */                      } /* end of component loop i */
1787                  } /* end of k0 loop */                  } /* end of k0 loop */
1788              } /* end of face 2 */              } /* end of face 2 */
1789              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1790  #pragma omp for nowait  #pragma omp for nowait
1791                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1792                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1793                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1794                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1795                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1796                          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];
1797                          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];
1798                      } /* end of component loop i */                      } /* end of component loop i */
1799                  } /* end of k0 loop */                  } /* end of k0 loop */
1800              } /* end of face 3 */              } /* end of face 3 */
1801          } // end of parallel section          } /* end of parallel section */
1802      }      }
1803  }  }
1804    
# Line 1466  void Rectangle::assemblePDESingle(Paso_S Line 1808  void Rectangle::assemblePDESingle(Paso_S
1808          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
1809          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
1810  {  {
1811      const double h0 = m_l0/m_gNE0;      /* GENERATOR SNIP_PDE_SINGLE_PRE TOP */
1812      const double h1 = m_l1/m_gNE1;      const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];
     const double w0 = -0.1555021169820365539*h1/h0;  
1813      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
1814      const double w2 = -0.15550211698203655390;      const double w2 = -0.15550211698203655390;
1815      const double w3 = 0.041666666666666666667*h0/h1;      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];
1816      const double w4 = 0.15550211698203655390;      const double w4 = -0.01116454968463011277*m_dx[1]/m_dx[0];
1817      const double w5 = -0.041666666666666666667;      const double w5 = 0.011164549684630112770;
1818      const double w6 = -0.01116454968463011277*h1/h0;      const double w6 = -0.041666666666666666667*m_dx[1]/m_dx[0];
1819      const double w7 = 0.011164549684630112770;      const double w7 = 0.1555021169820365539*m_dx[0]/m_dx[1];
1820      const double w8 = -0.011164549684630112770;      const double w8 = 0.01116454968463011277*m_dx[0]/m_dx[1];
1821      const double w9 = -0.041666666666666666667*h1/h0;      const double w9 = -1./4;
1822      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.16666666666666666667*m_dx[0]/m_dx[1];
1823      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = -0.032861463941450536761*m_dx[1];
1824      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = -0.032861463941450536761*m_dx[0];
1825      const double w13 = 0.01116454968463011277*h0/h1;      const double w13 = -0.12264065304058601714*m_dx[1];
1826      const double w14 = 0.01116454968463011277*h1/h0;      const double w14 = -0.0023593469594139828636*m_dx[1];
1827      const double w15 = 0.041666666666666666667*h1/h0;      const double w15 = -0.008805202725216129906*m_dx[0];
1828      const double w16 = -0.01116454968463011277*h0/h1;      const double w16 = -0.008805202725216129906*m_dx[1];
1829      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.12264065304058601714*m_dx[0];
1830      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.0023593469594139828636*m_dx[0];
1831      const double w19 = 0.25;      const double w19 = -0.16666666666666666667*m_dx[1];
1832      const double w20 = -0.25;      const double w20 = -0.083333333333333333333*m_dx[0];
1833      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.025917019497006092316*m_dx[0]*m_dx[1];
1834      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = 0.0018607582807716854616*m_dx[0]*m_dx[1];
1835      const double w23 = -0.16666666666666666667*h0/h1;      const double w23 = 0.0069444444444444444444*m_dx[0]*m_dx[1];
1836      const double w24 = 0.33333333333333333333*h1/h0;      const double w24 = 0.09672363354357992482*m_dx[0]*m_dx[1];
1837      const double w25 = 0.33333333333333333333*h0/h1;      const double w25 = 0.000498588678642297402*m_dx[0]*m_dx[1];
1838      const double w26 = 0.16666666666666666667*h1/h0;      const double w26 = 0.19716878364870322056*m_dx[1];
1839      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = 0.052831216351296779436*m_dx[1];
1840      const double w28 = -0.032861463941450536761*h1;      const double w28 = 0.19716878364870322056*m_dx[0];
1841      const double w29 = -0.032861463941450536761*h0;      const double w29 = 0.052831216351296779436*m_dx[0];
1842      const double w30 = -0.12264065304058601714*h1;      /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */
     const double w31 = -0.0023593469594139828636*h1;  
     const double w32 = -0.008805202725216129906*h0;  
     const double w33 = -0.008805202725216129906*h1;  
     const double w34 = 0.032861463941450536761*h1;  
     const double w35 = 0.008805202725216129906*h1;  
     const double w36 = 0.008805202725216129906*h0;  
     const double w37 = 0.0023593469594139828636*h1;  
     const double w38 = 0.12264065304058601714*h1;  
     const double w39 = 0.032861463941450536761*h0;  
     const double w40 = -0.12264065304058601714*h0;  
     const double w41 = -0.0023593469594139828636*h0;  
     const double w42 = 0.0023593469594139828636*h0;  
     const double w43 = 0.12264065304058601714*h0;  
     const double w44 = -0.16666666666666666667*h1;  
     const double w45 = -0.083333333333333333333*h0;  
     const double w46 = 0.083333333333333333333*h1;  
     const double w47 = 0.16666666666666666667*h1;  
     const double w48 = 0.083333333333333333333*h0;  
     const double w49 = -0.16666666666666666667*h0;  
     const double w50 = 0.16666666666666666667*h0;  
     const double w51 = -0.083333333333333333333*h1;  
     const double w52 = 0.025917019497006092316*h0*h1;  
     const double w53 = 0.0018607582807716854616*h0*h1;  
     const double w54 = 0.0069444444444444444444*h0*h1;  
     const double w55 = 0.09672363354357992482*h0*h1;  
     const double w56 = 0.00049858867864229740201*h0*h1;  
     const double w57 = 0.055555555555555555556*h0*h1;  
     const double w58 = 0.027777777777777777778*h0*h1;  
     const double w59 = 0.11111111111111111111*h0*h1;  
     const double w60 = -0.19716878364870322056*h1;  
     const double w61 = -0.19716878364870322056*h0;  
     const double w62 = -0.052831216351296779436*h0;  
     const double w63 = -0.052831216351296779436*h1;  
     const double w64 = 0.19716878364870322056*h1;  
     const double w65 = 0.052831216351296779436*h1;  
     const double w66 = 0.19716878364870322056*h0;  
     const double w67 = 0.052831216351296779436*h0;  
     const double w68 = -0.5*h1;  
     const double w69 = -0.5*h0;  
     const double w70 = 0.5*h1;  
     const double w71 = 0.5*h0;  
     const double w72 = 0.1555021169820365539*h0*h1;  
     const double w73 = 0.041666666666666666667*h0*h1;  
     const double w74 = 0.01116454968463011277*h0*h1;  
     const double w75 = 0.25*h0*h1;  
1843    
1844      rhs.requireWrite();      rhs.requireWrite();
1845  #pragma omp parallel  #pragma omp parallel
1846      {      {
1847          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1848  #pragma omp for  #pragma omp for
1849              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {              for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
1850                  for (index_t k0=0; k0<m_NE0; ++k0)  {                  for (index_t k0=0; k0<m_NE[0]; ++k0)  {
1851                      bool add_EM_S=false;                      bool add_EM_S=false;
1852                      bool add_EM_F=false;                      bool add_EM_F=false;
1853                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1854                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1855                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE[0]*k1;
1856                        /* GENERATOR SNIP_PDE_SINGLE TOP */
1857                      ///////////////                      ///////////////
1858                      // process A //                      // process A //
1859                      ///////////////                      ///////////////
1860                      if (!A.isEmpty()) {                      if (!A.isEmpty()) {
1861                          add_EM_S=true;                          add_EM_S = true;
1862                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p = const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1863                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
1864                              const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];                              const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
                             const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];  
1865                              const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];                              const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1866                                const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1867                              const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];                              const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1868                              const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];                              const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
                             const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];  
1869                              const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];                              const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1870                                const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1871                              const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];                              const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1872                              const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];                              const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
                             const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];  
1873                              const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];                              const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1874                                const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1875                              const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];                              const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1876                              const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];                              const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
                             const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];  
1877                              const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];                              const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1878                                const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1879                              const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];                              const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1880                              const double tmp0_0 = A_01_0 + A_01_3;                              const double tmp0 = w3*(A_11_0 + A_11_1 + A_11_2 + A_11_3);
1881                              const double tmp1_0 = A_00_0 + A_00_1;                              const double tmp1 = w1*(A_01_0 + A_01_3 - A_10_1 - A_10_2);
1882                              const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                              const double tmp2 = w4*(A_00_2 + A_00_3);
1883                              const double tmp3_0 = A_00_2 + A_00_3;                              const double tmp3 = w0*(A_00_0 + A_00_1);
1884                              const double tmp4_0 = A_10_1 + A_10_2;                              const double tmp4 = w5*(A_01_2 - A_10_3);
1885                              const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                              const double tmp5 = w2*(-A_01_1 + A_10_0);
1886                              const double tmp6_0 = A_01_3 + A_10_0;                              const double tmp6 = w5*(A_01_3 + A_10_0);
1887                              const double tmp7_0 = A_01_0 + A_10_3;                              const double tmp7 = w3*(-A_11_0 - A_11_1 - A_11_2 - A_11_3);
1888                              const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                              const double tmp8 = w6*(A_00_0 + A_00_1 + A_00_2 + A_00_3);
1889                              const double tmp9_0 = A_01_0 + A_10_0;                              const double tmp9 = w1*(A_01_1 + A_01_2 + A_10_1 + A_10_2);
1890                              const double tmp12_0 = A_11_0 + A_11_2;                              const double tmp10 = w2*(-A_01_0 - A_10_3);
1891                              const double tmp10_0 = A_01_3 + A_10_3;                              const double tmp11 = w4*(A_00_0 + A_00_1);
1892                              const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                              const double tmp12 = w0*(A_00_2 + A_00_3);
1893                              const double tmp13_0 = A_01_2 + A_10_1;                              const double tmp13 = w5*(A_01_1 - A_10_0);
1894                              const double tmp11_0 = A_11_1 + A_11_3;                              const double tmp14 = w2*(-A_01_2 + A_10_3);
1895                              const double tmp18_0 = A_01_1 + A_10_1;                              const double tmp15 = w7*(A_11_0 + A_11_2);
1896                              const double tmp15_0 = A_01_1 + A_10_2;                              const double tmp16 = w4*(-A_00_2 - A_00_3);
1897                              const double tmp16_0 = A_10_0 + A_10_3;                              const double tmp17 = w0*(-A_00_0 - A_00_1);
1898                              const double tmp17_0 = A_01_1 + A_01_2;                              const double tmp18 = w5*(A_01_3 + A_10_3);
1899                              const double tmp19_0 = A_01_2 + A_10_2;                              const double tmp19 = w8*(A_11_1 + A_11_3);
1900                              const double tmp0_1 = A_10_3*w8;                              const double tmp20 = w2*(-A_01_0 - A_10_0);
1901                              const double tmp1_1 = tmp0_0*w1;                              const double tmp21 = w7*(A_11_1 + A_11_3);
1902                              const double tmp2_1 = A_01_1*w4;                              const double tmp22 = w4*(-A_00_0 - A_00_1);
1903                              const double tmp3_1 = tmp1_0*w0;                              const double tmp23 = w0*(-A_00_2 - A_00_3);
1904                              const double tmp4_1 = A_01_2*w7;                              const double tmp24 = w5*(A_01_0 + A_10_0);
1905                              const double tmp5_1 = tmp2_0*w3;                              const double tmp25 = w8*(A_11_0 + A_11_2);
1906                              const double tmp6_1 = tmp3_0*w6;                              const double tmp26 = w2*(-A_01_3 - A_10_3);
1907                              const double tmp7_1 = A_10_0*w2;                              const double tmp27 = w5*(-A_01_1 - A_10_2);
1908                              const double tmp8_1 = tmp4_0*w5;                              const double tmp28 = w1*(-A_01_0 - A_01_3 - A_10_0 - A_10_3);
1909                              const double tmp9_1 = tmp2_0*w10;                              const double tmp29 = w2*(A_01_2 + A_10_1);
1910                              const double tmp14_1 = A_10_0*w8;                              const double tmp30 = w7*(-A_11_1 - A_11_3);
1911                              const double tmp23_1 = tmp3_0*w14;                              const double tmp31 = w1*(-A_01_1 - A_01_2 + A_10_0 + A_10_3);
1912                              const double tmp35_1 = A_01_0*w8;                              const double tmp32 = w5*(-A_01_0 + A_10_2);
1913                              const double tmp54_1 = tmp13_0*w8;                              const double tmp33 = w8*(-A_11_0 - A_11_2);
1914                              const double tmp20_1 = tmp9_0*w4;                              const double tmp34 = w6*(-A_00_0 - A_00_1 - A_00_2 - A_00_3);
1915                              const double tmp25_1 = tmp12_0*w12;                              const double tmp35 = w2*(A_01_3 - A_10_1);
1916                              const double tmp44_1 = tmp7_0*w7;                              const double tmp36 = w5*(A_01_0 + A_10_3);
1917                              const double tmp26_1 = tmp10_0*w4;                              const double tmp37 = w2*(-A_01_3 - A_10_0);
1918                              const double tmp52_1 = tmp18_0*w8;                              const double tmp38 = w7*(-A_11_0 - A_11_2);
1919                              const double tmp48_1 = A_10_1*w7;                              const double tmp39 = w5*(-A_01_3 + A_10_1);
1920                              const double tmp46_1 = A_01_3*w8;                              const double tmp40 = w8*(-A_11_1 - A_11_3);
1921                              const double tmp50_1 = A_01_0*w2;                              const double tmp41 = w2*(A_01_0 - A_10_2);
1922                              const double tmp56_1 = tmp19_0*w8;                              const double tmp42 = w5*(A_01_1 - A_10_3);
1923                              const double tmp19_1 = A_10_3*w2;                              const double tmp43 = w2*(-A_01_2 + A_10_0);
1924                              const double tmp47_1 = A_10_2*w4;                              const double tmp44 = w5*(A_01_2 - A_10_0);
1925                              const double tmp16_1 = tmp3_0*w0;                              const double tmp45 = w2*(-A_01_1 + A_10_3);
1926                              const double tmp18_1 = tmp1_0*w6;                              const double tmp46 = w5*(-A_01_0 + A_10_1);
1927                              const double tmp31_1 = tmp11_0*w12;                              const double tmp47 = w2*(A_01_3 - A_10_2);
1928                              const double tmp55_1 = tmp15_0*w2;                              const double tmp48 = w5*(-A_01_1 - A_10_1);
1929                              const double tmp39_1 = A_10_2*w7;                              const double tmp49 = w2*(A_01_2 + A_10_2);
1930                              const double tmp11_1 = tmp6_0*w7;                              const double tmp50 = w5*(-A_01_3 + A_10_2);
1931                              const double tmp40_1 = tmp11_0*w17;                              const double tmp51 = w2*(A_01_0 - A_10_1);
1932                              const double tmp34_1 = tmp15_0*w8;                              const double tmp52 = w5*(-A_01_2 - A_10_1);
1933                              const double tmp33_1 = tmp14_0*w5;                              const double tmp53 = w2*(A_01_1 + A_10_2);
1934                              const double tmp24_1 = tmp11_0*w13;                              const double tmp54 = w5*(-A_01_2 - A_10_2);
1935                              const double tmp43_1 = tmp17_0*w5;                              const double tmp55 = w2*(A_01_1 + A_10_1);
1936                              const double tmp15_1 = A_01_2*w4;                              EM_S[INDEX2(0,0,4)]+=tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp9;
1937                              const double tmp53_1 = tmp19_0*w2;                              EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;
1938                              const double tmp27_1 = tmp3_0*w11;                              EM_S[INDEX2(0,2,4)]+=tmp31 + tmp34 + tmp38 + tmp39 + tmp40 + tmp41;
1939                              const double tmp32_1 = tmp13_0*w2;                              EM_S[INDEX2(0,3,4)]+=tmp28 + tmp52 + tmp53 + tmp7 + tmp8;
1940                              const double tmp10_1 = tmp5_0*w9;                              EM_S[INDEX2(1,0,4)]+=tmp0 + tmp2 + tmp3 + tmp31 + tmp50 + tmp51;
1941                              const double tmp37_1 = A_10_1*w4;                              EM_S[INDEX2(1,1,4)]+=tmp16 + tmp17 + tmp21 + tmp25 + tmp28 + tmp54 + tmp55;
1942                              const double tmp38_1 = tmp5_0*w15;                              EM_S[INDEX2(1,2,4)]+=tmp10 + tmp6 + tmp7 + tmp8 + tmp9;
1943                              const double tmp17_1 = A_01_1*w7;                              EM_S[INDEX2(1,3,4)]+=tmp1 + tmp30 + tmp33 + tmp34 + tmp44 + tmp45;
1944                              const double tmp12_1 = tmp7_0*w4;                              EM_S[INDEX2(2,0,4)]+=tmp1 + tmp34 + tmp38 + tmp40 + tmp42 + tmp43;
1945                              const double tmp22_1 = tmp10_0*w7;                              EM_S[INDEX2(2,1,4)]+=tmp36 + tmp37 + tmp7 + tmp8 + tmp9;
1946                              const double tmp57_1 = tmp18_0*w2;                              EM_S[INDEX2(2,2,4)]+=tmp15 + tmp19 + tmp22 + tmp23 + tmp28 + tmp48 + tmp49;
1947                              const double tmp28_1 = tmp9_0*w7;                              EM_S[INDEX2(2,3,4)]+=tmp0 + tmp11 + tmp12 + tmp31 + tmp46 + tmp47;
1948                              const double tmp29_1 = tmp1_0*w14;                              EM_S[INDEX2(3,0,4)]+=tmp27 + tmp28 + tmp29 + tmp7 + tmp8;
1949                              const double tmp51_1 = tmp11_0*w16;                              EM_S[INDEX2(3,1,4)]+=tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35;
1950                              const double tmp42_1 = tmp12_0*w16;                              EM_S[INDEX2(3,2,4)]+=tmp0 + tmp1 + tmp11 + tmp12 + tmp13 + tmp14;
1951                              const double tmp49_1 = tmp12_0*w17;                              EM_S[INDEX2(3,3,4)]+=tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp9;
                             const double tmp21_1 = tmp1_0*w11;  
                             const double tmp45_1 = tmp6_0*w4;  
                             const double tmp13_1 = tmp8_0*w1;  
                             const double tmp36_1 = tmp16_0*w1;  
                             const double tmp41_1 = A_01_3*w2;  
                             const double tmp30_1 = tmp12_0*w13;  
                             EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
1952                          } else { // constant data                          } else { // constant data
1953                              const double A_00 = A_p[INDEX2(0,0,2)];                              const double A_00 = A_p[INDEX2(0,0,2)];
                             const double A_10 = A_p[INDEX2(1,0,2)];  
1954                              const double A_01 = A_p[INDEX2(0,1,2)];                              const double A_01 = A_p[INDEX2(0,1,2)];
1955                                const double A_10 = A_p[INDEX2(1,0,2)];
1956                              const double A_11 = A_p[INDEX2(1,1,2)];                              const double A_11 = A_p[INDEX2(1,1,2)];
1957                              const double tmp0_0 = A_01 + A_10;                              const double tmp0 = w9*(-A_01 - A_10);
1958                              const double tmp0_1 = A_00*w18;                              const double tmp1 = w9*(-A_01 + A_10);
1959                              const double tmp1_1 = A_01*w19;                              const double tmp2 = w9*(A_01 + A_10);
1960                              const double tmp2_1 = A_10*w20;                              const double tmp3 = w9*(A_01 - A_10);
1961                              const double tmp3_1 = A_11*w21;                              EM_S[INDEX2(0,0,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;
1962                              const double tmp4_1 = A_00*w22;                              EM_S[INDEX2(0,1,4)]+=8*A_00*w6 + 6*A_01*w1 + A_10*w9 + 4*A_11*w3;
1963                              const double tmp5_1 = tmp0_0*w19;                              EM_S[INDEX2(0,2,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;
1964                              const double tmp6_1 = A_11*w23;                              EM_S[INDEX2(0,3,4)]+=4*A_00*w6 + A_11*w10 + tmp2;
1965                              const double tmp7_1 = A_11*w25;                              EM_S[INDEX2(1,0,4)]+=8*A_00*w6 - A_11*w10 + tmp3;
1966                              const double tmp8_1 = A_00*w24;                              EM_S[INDEX2(1,1,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;
1967                              const double tmp9_1 = tmp0_0*w20;                              EM_S[INDEX2(1,2,4)]+=4*A_00*w6 + A_11*w10 + tmp0;
1968                              const double tmp10_1 = A_01*w20;                              EM_S[INDEX2(1,3,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;
1969                              const double tmp11_1 = A_11*w27;                              EM_S[INDEX2(2,0,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;
1970                              const double tmp12_1 = A_00*w26;                              EM_S[INDEX2(2,1,4)]+=4*A_00*w6 + A_11*w10 + tmp0;
1971                              const double tmp13_1 = A_10*w19;                              EM_S[INDEX2(2,2,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;
1972                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(2,3,4)]+=8*A_00*w6 - A_11*w10 + tmp3;
1973                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                              EM_S[INDEX2(3,0,4)]+=4*A_00*w6 + A_11*w10 + tmp2;
1974                              EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(3,1,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;
1975                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(3,2,4)]+=8*A_00*w6 - A_11*w10 + tmp1;
1976                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(3,3,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;
                             EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
1977                          }                          }
1978                      }                      }
1979                      ///////////////                      ///////////////
# Line 1727  void Rectangle::assemblePDESingle(Paso_S Line 1991  void Rectangle::assemblePDESingle(Paso_S
1991                              const double B_1_2 = B_p[INDEX2(1,2,2)];                              const double B_1_2 = B_p[INDEX2(1,2,2)];
1992                              const double B_0_3 = B_p[INDEX2(0,3,2)];                              const double B_0_3 = B_p[INDEX2(0,3,2)];
1993                              const double B_1_3 = B_p[INDEX2(1,3,2)];                              const double B_1_3 = B_p[INDEX2(1,3,2)];
1994                              const double tmp0_0 = B_1_0 + B_1_1;                              const double tmp0 = w15*(B_1_2 + B_1_3);
1995                              const double tmp1_0 = B_1_2 + B_1_3;                              const double tmp1 = w12*(B_1_0 + B_1_1);
1996                              const double tmp2_0 = B_0_1 + B_0_3;                              const double tmp2 = w15*(B_1_0 + B_1_1);
1997                              const double tmp3_0 = B_0_0 + B_0_2;                              const double tmp3 = w16*(-B_0_1 - B_0_3);
1998                              const double tmp63_1 = B_1_1*w42;                              const double tmp4 = w11*(-B_0_0 - B_0_2);
1999                              const double tmp79_1 = B_1_1*w40;                              const double tmp5 = w12*(B_1_2 + B_1_3);
2000                              const double tmp37_1 = tmp3_0*w35;                              const double tmp6 = w15*(-B_1_0 - B_1_1);
2001                              const double tmp8_1 = tmp0_0*w32;                              const double tmp7 = w12*(-B_1_2 - B_1_3);
2002                              const double tmp71_1 = B_0_1*w34;                              const double tmp8 = w15*(-B_1_2 - B_1_3);
2003                              const double tmp19_1 = B_0_3*w31;                              const double tmp9 = w12*(-B_1_0 - B_1_1);
2004                              const double tmp15_1 = B_0_3*w34;                              const double tmp10 = w11*(-B_0_1 - B_0_3);
2005                              const double tmp9_1 = tmp3_0*w34;                              const double tmp11 = w16*(-B_0_0 - B_0_2);
2006                              const double tmp35_1 = B_1_0*w36;                              const double tmp12 = w16*(B_0_0 + B_0_2);
2007                              const double tmp66_1 = B_0_3*w28;                              const double tmp13 = w11*(B_0_1 + B_0_3);
2008                              const double tmp28_1 = B_1_0*w42;                              const double tmp14 = w11*(B_0_0 + B_0_2);
2009                              const double tmp22_1 = B_1_0*w40;                              const double tmp15 = w16*(B_0_1 + B_0_3);
2010                              const double tmp16_1 = B_1_2*w29;                              EM_S[INDEX2(0,0,4)]+=B_0_0*w13 + B_0_1*w11 + B_0_2*w16 + B_0_3*w14 + B_1_0*w17 + B_1_1*w15 + B_1_2*w12 + B_1_3*w18;
2011                              const double tmp6_1 = tmp2_0*w35;                              EM_S[INDEX2(0,1,4)]+=B_0_0*w11 + B_0_1*w13 + B_0_2*w14 + B_0_3*w16 + tmp0 + tmp1;
2012                              const double tmp55_1 = B_1_3*w40;                              EM_S[INDEX2(0,2,4)]+=B_1_0*w12 + B_1_1*w18 + B_1_2*w17 + B_1_3*w15 + tmp14 + tmp15;
2013                              const double tmp50_1 = B_1_3*w42;                              EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp2 + tmp5;
2014                              const double tmp7_1 = tmp1_0*w29;                              EM_S[INDEX2(1,0,4)]+=-B_0_0*w13 - B_0_1*w11 - B_0_2*w16 - B_0_3*w14 + tmp0 + tmp1;
2015                              const double tmp1_1 = tmp1_0*w32;                              EM_S[INDEX2(1,1,4)]+=-B_0_0*w11 - B_0_1*w13 - B_0_2*w14 - B_0_3*w16 + B_1_0*w15 + B_1_1*w17 + B_1_2*w18 + B_1_3*w12;
2016                              const double tmp57_1 = B_0_3*w30;                              EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2017                              const double tmp18_1 = B_1_1*w32;                              EM_S[INDEX2(1,3,4)]+=B_1_0*w18 + B_1_1*w12 + B_1_2*w15 + B_1_3*w17 + tmp10 + tmp11;
2018                              const double tmp53_1 = B_1_0*w41;                              EM_S[INDEX2(2,0,4)]+=-B_1_0*w17 - B_1_1*w15 - B_1_2*w12 - B_1_3*w18 + tmp14 + tmp15;
2019                              const double tmp61_1 = B_1_3*w36;                              EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2020                              const double tmp27_1 = B_0_3*w38;                              EM_S[INDEX2(2,2,4)]+=B_0_0*w16 + B_0_1*w14 + B_0_2*w13 + B_0_3*w11 - B_1_0*w12 - B_1_1*w18 - B_1_2*w17 - B_1_3*w15;
2021                              const double tmp64_1 = B_0_2*w30;                              EM_S[INDEX2(2,3,4)]+=B_0_0*w14 + B_0_1*w16 + B_0_2*w11 + B_0_3*w13 + tmp6 + tmp7;
2022                              const double tmp76_1 = B_0_1*w38;                              EM_S[INDEX2(3,0,4)]+=tmp3 + tmp4 + tmp8 + tmp9;
2023                              const double tmp39_1 = tmp2_0*w34;                              EM_S[INDEX2(3,1,4)]+=-B_1_0*w15 - B_1_1*w17 - B_1_2*w18 - B_1_3*w12 + tmp10 + tmp11;
2024                              const double tmp62_1 = B_0_1*w31;                              EM_S[INDEX2(3,2,4)]+=-B_0_0*w16 - B_0_1*w14 - B_0_2*w13 - B_0_3*w11 + tmp6 + tmp7;
2025                              const double tmp56_1 = B_0_0*w31;                              EM_S[INDEX2(3,3,4)]+=-B_0_0*w14 - B_0_1*w16 - B_0_2*w11 - B_0_3*w13 - B_1_0*w18 - B_1_1*w12 - B_1_2*w15 - B_1_3*w17;
                             const double tmp49_1 = B_1_1*w36;  
                             const double tmp2_1 = B_0_2*w31;  
                             const double tmp23_1 = B_0_2*w33;  
                             const double tmp38_1 = B_1_1*w43;  
                             const double tmp74_1 = B_1_2*w41;  
                             const double tmp43_1 = B_1_1*w41;  
                             const double tmp58_1 = B_0_2*w28;  
                             const double tmp67_1 = B_0_0*w33;  
                             const double tmp33_1 = tmp0_0*w39;  
                             const double tmp4_1 = B_0_0*w28;  
                             const double tmp20_1 = B_0_0*w30;  
                             const double tmp13_1 = B_0_2*w38;  
                             const double tmp65_1 = B_1_2*w43;  
                             const double tmp0_1 = tmp0_0*w29;  
                             const double tmp41_1 = tmp3_0*w33;  
                             const double tmp73_1 = B_0_2*w37;  
                             const double tmp69_1 = B_0_0*w38;  
                             const double tmp48_1 = B_1_2*w39;  
                             const double tmp59_1 = B_0_1*w33;  
                             const double tmp17_1 = B_1_3*w41;  
                             const double tmp5_1 = B_0_3*w33;  
                             const double tmp3_1 = B_0_1*w30;  
                             const double tmp21_1 = B_0_1*w28;  
                             const double tmp42_1 = B_1_0*w29;  
                             const double tmp54_1 = B_1_2*w32;  
                             const double tmp60_1 = B_1_0*w39;  
                             const double tmp32_1 = tmp1_0*w36;  
                             const double tmp10_1 = B_0_1*w37;  
                             const double tmp14_1 = B_0_0*w35;  
                             const double tmp29_1 = B_0_1*w35;  
                             const double tmp26_1 = B_1_2*w36;  
                             const double tmp30_1 = B_1_3*w43;  
                             const double tmp70_1 = B_0_2*w35;  
                             const double tmp34_1 = B_1_3*w39;  
                             const double tmp51_1 = B_1_0*w43;  
                             const double tmp31_1 = B_0_2*w34;  
                             const double tmp45_1 = tmp3_0*w28;  
                             const double tmp11_1 = tmp1_0*w39;  
                             const double tmp52_1 = B_1_1*w29;  
                             const double tmp44_1 = B_1_3*w32;  
                             const double tmp25_1 = B_1_1*w39;  
                             const double tmp47_1 = tmp2_0*w33;  
                             const double tmp72_1 = B_1_3*w29;  
                             const double tmp40_1 = tmp2_0*w28;  
                             const double tmp46_1 = B_1_2*w40;  
                             const double tmp36_1 = B_1_2*w42;  
                             const double tmp24_1 = B_0_0*w37;  
                             const double tmp77_1 = B_0_3*w35;  
                             const double tmp68_1 = B_0_3*w37;  
                             const double tmp78_1 = B_0_0*w34;  
                             const double tmp12_1 = tmp0_0*w36;  
                             const double tmp75_1 = B_1_0*w32;  
                             EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
2026                          } else { // constant data                          } else { // constant data
2027                              const double B_0 = B_p[0];                              const double B_0 = B_p[0];
2028                              const double B_1 = B_p[1];                              const double B_1 = B_p[1];
2029                              const double tmp0_1 = B_0*w44;                              EM_S[INDEX2(0,0,4)]+=B_0*w19 + 2*B_1*w20;
2030                              const double tmp1_1 = B_1*w45;                              EM_S[INDEX2(0,1,4)]+=B_0*w19 + B_1*w20;
2031                              const double tmp2_1 = B_0*w46;                              EM_S[INDEX2(0,2,4)]+=B_0*w19/2 + 2*B_1*w20;
2032                              const double tmp3_1 = B_0*w47;                              EM_S[INDEX2(0,3,4)]+=B_0*w19/2 + B_1*w20;
2033                              const double tmp4_1 = B_1*w48;                              EM_S[INDEX2(1,0,4)]+=-B_0*w19 + B_1*w20;
2034                              const double tmp5_1 = B_1*w49;                              EM_S[INDEX2(1,1,4)]+=-B_0*w19 + 2*B_1*w20;
2035                              const double tmp6_1 = B_1*w50;                              EM_S[INDEX2(1,2,4)]+=-B_0*w19/2 + B_1*w20;
2036                              const double tmp7_1 = B_0*w51;                              EM_S[INDEX2(1,3,4)]+=-B_0*w19/2 + 2*B_1*w20;
2037                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(2,0,4)]+=B_0*w19/2 - 2*B_1*w20;
2038                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(2,1,4)]+=B_0*w19/2 - B_1*w20;
2039                              EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;                              EM_S[INDEX2(2,2,4)]+=B_0*w19 - 2*B_1*w20;
2040                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;                              EM_S[INDEX2(2,3,4)]+=B_0*w19 - B_1*w20;
2041                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(3,0,4)]+=-B_0*w19/2 - B_1*w20;
2042                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(3,1,4)]+=-B_0*w19/2 - 2*B_1*w20;
2043                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;                              EM_S[INDEX2(3,2,4)]+=-B_0*w19 - B_1*w20;
2044                              EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;                              EM_S[INDEX2(3,3,4)]+=-B_0*w19 - 2*B_1*w20;
                             EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;  
2045                          }                          }
2046                      }                      }
2047                      ///////////////                      ///////////////
# Line 1871  void Rectangle::assemblePDESingle(Paso_S Line 2059  void Rectangle::assemblePDESingle(Paso_S
2059                              const double C_1_2 = C_p[INDEX2(1,2,2)];                              const double C_1_2 = C_p[INDEX2(1,2,2)];
2060                              const double C_0_3 = C_p[INDEX2(0,3,2)];                              const double C_0_3 = C_p[INDEX2(0,3,2)];
2061                              const double C_1_3 = C_p[INDEX2(1,3,2)];                              const double C_1_3 = C_p[INDEX2(1,3,2)];
2062                              const double tmp0_0 = C_1_0 + C_1_1;                              const double tmp0 = w15*(C_1_2 + C_1_3);
2063                              const double tmp1_0 = C_1_2 + C_1_3;                              const double tmp1 = w12*(C_1_0 + C_1_1);
2064                              const double tmp2_0 = C_0_1 + C_0_3;                              const double tmp2 = w15*(-C_1_2 - C_1_3);
2065                              const double tmp3_0 = C_0_0 + C_0_2;                              const double tmp3 = w16*(C_0_0 + C_0_2);
2066                              const double tmp64_1 = C_0_2*w30;                              const double tmp4 = w11*(C_0_1 + C_0_3);
2067                              const double tmp14_1 = C_0_2*w28;                              const double tmp5 = w12*(-C_1_0 - C_1_1);
2068                              const double tmp19_1 = C_0_3*w31;                              const double tmp6 = w15*(-C_1_0 - C_1_1);
2069                              const double tmp22_1 = C_1_0*w40;                              const double tmp7 = w12*(-C_1_2 - C_1_3);
2070                              const double tmp37_1 = tmp3_0*w35;                              const double tmp8 = w15*(C_1_0 + C_1_1);
2071                              const double tmp29_1 = C_0_1*w35;                              const double tmp9 = w12*(C_1_2 + C_1_3);
2072                              const double tmp73_1 = C_0_2*w37;                              const double tmp10 = w11*(-C_0_1 - C_0_3);
2073                              const double tmp74_1 = C_1_2*w41;                              const double tmp11 = w16*(-C_0_0 - C_0_2);
2074                              const double tmp52_1 = C_1_3*w39;                              const double tmp12 = w16*(-C_0_1 - C_0_3);
2075                              const double tmp25_1 = C_1_1*w39;                              const double tmp13 = w11*(-C_0_0 - C_0_2);
2076                              const double tmp62_1 = C_0_1*w31;                              const double tmp14 = w11*(C_0_0 + C_0_2);
2077                              const double tmp79_1 = C_1_1*w40;                              const double tmp15 = w16*(C_0_1 + C_0_3);
2078                              const double tmp43_1 = C_1_1*w36;                              EM_S[INDEX2(0,0,4)]+=C_0_0*w13 + C_0_1*w11 + C_0_2*w16 + C_0_3*w14 + C_1_0*w17 + C_1_1*w15 + C_1_2*w12 + C_1_3*w18;
2079                              const double tmp27_1 = C_0_3*w38;                              EM_S[INDEX2(0,1,4)]+=-C_0_0*w13 - C_0_1*w11 - C_0_2*w16 - C_0_3*w14 + tmp0 + tmp1;
2080                              const double tmp28_1 = C_1_0*w42;                              EM_S[INDEX2(0,2,4)]+=-C_1_0*w17 - C_1_1*w15 - C_1_2*w12 - C_1_3*w18 + tmp14 + tmp15;
2081                              const double tmp63_1 = C_1_1*w42;                              EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp2 + tmp5;
2082                              const double tmp59_1 = C_0_3*w34;                              EM_S[INDEX2(1,0,4)]+=C_0_0*w11 + C_0_1*w13 + C_0_2*w14 + C_0_3*w16 + tmp0 + tmp1;
2083                              const double tmp72_1 = C_1_3*w29;                              EM_S[INDEX2(1,1,4)]+=-C_0_0*w11 - C_0_1*w13 - C_0_2*w14 - C_0_3*w16 + C_1_0*w15 + C_1_1*w17 + C_1_2*w18 + C_1_3*w12;
2084                              const double tmp40_1 = tmp2_0*w35;                              EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2085                              const double tmp13_1 = C_0_3*w30;                              EM_S[INDEX2(1,3,4)]+=-C_1_0*w15 - C_1_1*w17 - C_1_2*w18 - C_1_3*w12 + tmp10 + tmp11;
2086                              const double tmp51_1 = C_1_2*w40;                              EM_S[INDEX2(2,0,4)]+=C_1_0*w12 + C_1_1*w18 + C_1_2*w17 + C_1_3*w15 + tmp14 + tmp15;
2087                              const double tmp54_1 = C_1_2*w42;                              EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2088                              const double tmp12_1 = C_0_0*w31;                              EM_S[INDEX2(2,2,4)]+=C_0_0*w16 + C_0_1*w14 + C_0_2*w13 + C_0_3*w11 - C_1_0*w12 - C_1_1*w18 - C_1_2*w17 - C_1_3*w15;
2089                              const double tmp2_1 = tmp1_0*w32;                              EM_S[INDEX2(2,3,4)]+=-C_0_0*w16 - C_0_1*w14 - C_0_2*w13 - C_0_3*w11 + tmp6 + tmp7;
2090                              const double tmp68_1 = C_0_2*w31;                              EM_S[INDEX2(3,0,4)]+=tmp3 + tmp4 + tmp8 + tmp9;
2091                              const double tmp75_1 = C_1_0*w32;                              EM_S[INDEX2(3,1,4)]+=C_1_0*w18 + C_1_1*w12 + C_1_2*w15 + C_1_3*w17 + tmp10 + tmp11;
2092                              const double tmp49_1 = C_1_1*w41;                              EM_S[INDEX2(3,2,4)]+=C_0_0*w14 + C_0_1*w16 + C_0_2*w11 + C_0_3*w13 + tmp6 + tmp7;
2093                              const double tmp4_1 = C_0_2*w35;                              EM_S[INDEX2(3,3,4)]+=-C_0_0*w14 - C_0_1*w16 - C_0_2*w11 - C_0_3*w13 - C_1_0*w18 - C_1_1*w12 - C_1_2*w15 - C_1_3*w17;
                             const double tmp66_1 = C_0_3*w28;  
                             const double tmp56_1 = C_0_1*w37;  
                             const double tmp5_1 = C_0_1*w34;  
                             const double tmp38_1 = tmp2_0*w34;  
                             const double tmp76_1 = C_0_1*w38;  
                             const double tmp21_1 = C_0_1*w28;  
                             const double tmp69_1 = C_0_1*w30;  
                             const double tmp53_1 = C_1_0*w36;  
                             const double tmp42_1 = C_1_2*w39;  
                             const double tmp32_1 = tmp1_0*w29;  
                             const double tmp45_1 = C_1_0*w43;  
                             const double tmp33_1 = tmp0_0*w32;  
                             const double tmp35_1 = C_1_0*w41;  
                             const double tmp26_1 = C_1_2*w36;  
                             const double tmp67_1 = C_0_0*w33;  
                             const double tmp31_1 = C_0_2*w34;  
                             const double tmp20_1 = C_0_0*w30;  
                             const double tmp70_1 = C_0_0*w28;  
                             const double tmp8_1 = tmp0_0*w39;  
                             const double tmp30_1 = C_1_3*w43;  
                             const double tmp0_1 = tmp0_0*w29;  
                             const double tmp17_1 = C_1_3*w41;  
                             const double tmp58_1 = C_0_0*w35;  
                             const double tmp9_1 = tmp3_0*w33;  
                             const double tmp61_1 = C_1_3*w36;  
                             const double tmp41_1 = tmp3_0*w34;  
                             const double tmp50_1 = C_1_3*w32;  
                             const double tmp18_1 = C_1_1*w32;  
                             const double tmp6_1 = tmp1_0*w36;  
                             const double tmp3_1 = C_0_0*w38;  
                             const double tmp34_1 = C_1_1*w29;  
                             const double tmp77_1 = C_0_3*w35;  
                             const double tmp65_1 = C_1_2*w43;  
                             const double tmp71_1 = C_0_3*w33;  
                             const double tmp55_1 = C_1_1*w43;  
                             const double tmp46_1 = tmp3_0*w28;  
                             const double tmp24_1 = C_0_0*w37;  
                             const double tmp10_1 = tmp1_0*w39;  
                             const double tmp48_1 = C_1_0*w29;  
                             const double tmp15_1 = C_0_1*w33;  
                             const double tmp36_1 = C_1_2*w32;  
                             const double tmp60_1 = C_1_0*w39;  
                             const double tmp47_1 = tmp2_0*w33;  
                             const double tmp16_1 = C_1_2*w29;  
                             const double tmp1_1 = C_0_3*w37;  
                             const double tmp7_1 = tmp2_0*w28;  
                             const double tmp39_1 = C_1_3*w40;  
                             const double tmp44_1 = C_1_3*w42;  
                             const double tmp57_1 = C_0_2*w38;  
                             const double tmp78_1 = C_0_0*w34;  
                             const double tmp11_1 = tmp0_0*w36;  
                             const double tmp23_1 = C_0_2*w33;  
                             EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
2094                          } else { // constant data                          } else { // constant data
2095                              const double C_0 = C_p[0];                              const double C_0 = C_p[0];
2096                              const double C_1 = C_p[1];                              const double C_1 = C_p[1];
2097                              const double tmp0_1 = C_0*w47;                              EM_S[INDEX2(0,0,4)]+=C_0*w19 + 2*C_1*w20;
2098                              const double tmp1_1 = C_1*w45;                              EM_S[INDEX2(0,1,4)]+=-C_0*w19 + C_1*w20;
2099                              const double tmp2_1 = C_1*w48;                              EM_S[INDEX2(0,2,4)]+=C_0*w19/2 - 2*C_1*w20;
2100                              const double tmp3_1 = C_0*w51;                              EM_S[INDEX2(0,3,4)]+=-C_0*w19/2 - C_1*w20;
2101                              const double tmp4_1 = C_0*w44;                              EM_S[INDEX2(1,0,4)]+=C_0*w19 + C_1*w20;
2102                              const double tmp5_1 = C_1*w49;                              EM_S[INDEX2(1,1,4)]+=-C_0*w19 + 2*C_1*w20;
2103                              const double tmp6_1 = C_1*w50;                              EM_S[INDEX2(1,2,4)]+=C_0*w19/2 - C_1*w20;
2104                              const double tmp7_1 = C_0*w46;                              EM_S[INDEX2(1,3,4)]+=-C_0*w19/2 - 2*C_1*w20;
2105                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;                              EM_S[INDEX2(2,0,4)]+=C_0*w19/2 + 2*C_1*w20;
2106                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;                              EM_S[INDEX2(2,1,4)]+=-C_0*w19/2 + C_1*w20;
2107                              EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(2,2,4)]+=C_0*w19 - 2*C_1*w20;
2108                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(2,3,4)]+=-C_0*w19 - C_1*w20;
2109                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(3,0,4)]+=C_0*w19/2 + C_1*w20;
2110                              EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(3,1,4)]+=-C_0*w19/2 + 2*C_1*w20;
2111                              EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;                              EM_S[INDEX2(3,2,4)]+=C_0*w19 - C_1*w20;
2112                              EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;                              EM_S[INDEX2(3,3,4)]+=-C_0*w19 - 2*C_1*w20;
                             EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;  
2113                          }                          }
2114                      }                      }
2115                      ///////////////                      ///////////////
# Line 2011  void Rectangle::assemblePDESingle(Paso_S Line 2123  void Rectangle::assemblePDESingle(Paso_S
2123                              const double D_1 = D_p[1];                              const double D_1 = D_p[1];
2124                              const double D_2 = D_p[2];                              const double D_2 = D_p[2];
2125                              const double D_3 = D_p[3];                              const double D_3 = D_p[3];
2126                              const double tmp0_0 = D_0 + D_1;                              const double tmp0 = w21*(D_0 + D_1);
2127                              const double tmp1_0 = D_2 + D_3;                              const double tmp1 = w22*(D_2 + D_3);
2128                              const double tmp2_0 = D_0 + D_1 + D_2 + D_3;                              const double tmp2 = w23*(D_0 + D_1 + D_2 + D_3);
2129                              const double tmp3_0 = D_1 + D_2;                              const double tmp3 = w21*(D_2 + D_3);
2130                              const double tmp4_0 = D_1 + D_3;                              const double tmp4 = w22*(D_0 + D_1);
2131                              const double tmp5_0 = D_0 + D_2;                              const double tmp5 = w23*(D_1 + D_2);
2132                              const double tmp6_0 = D_0 + D_3;                              const double tmp6 = w21*(D_1 + D_3);
2133                              const double tmp0_1 = tmp0_0*w52;                              const double tmp7 = w22*(D_0 + D_2);
2134                              const double tmp1_1 = tmp1_0*w53;                              const double tmp8 = w21*(D_0 + D_2);
2135                              const double tmp2_1 = tmp2_0*w54;                              const double tmp9 = w22*(D_1 + D_3);
2136                              const double tmp3_1 = tmp1_0*w52;                              const double tmp10 = w23*(D_0 + D_3);
2137                              const double tmp4_1 = tmp0_0*w53;                              EM_S[INDEX2(0,0,4)]+=D_0*w24 + D_3*w25 + tmp5;
2138                              const double tmp5_1 = tmp3_0*w54;                              EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1;
2139                              const double tmp6_1 = D_0*w55;                              EM_S[INDEX2(0,2,4)]+=tmp8 + tmp9;
2140                              const double tmp7_1 = D_3*w56;                              EM_S[INDEX2(0,3,4)]+=tmp2;
2141                              const double tmp8_1 = D_3*w55;                              EM_S[INDEX2(1,0,4)]+=tmp0 + tmp1;
2142                              const double tmp9_1 = D_0*w56;                              EM_S[INDEX2(1,1,4)]+=D_1*w24 + D_2*w25 + tmp10;
2143                              const double tmp10_1 = tmp4_0*w52;                              EM_S[INDEX2(1,2,4)]+=tmp2;
2144                              const double tmp11_1 = tmp5_0*w53;                              EM_S[INDEX2(1,3,4)]+=tmp6 + tmp7;
2145                              const double tmp12_1 = tmp5_0*w52;                              EM_S[INDEX2(2,0,4)]+=tmp8 + tmp9;
2146                              const double tmp13_1 = tmp4_0*w53;                              EM_S[INDEX2(2,1,4)]+=tmp2;
2147                              const double tmp14_1 = tmp6_0*w54;                              EM_S[INDEX2(2,2,4)]+=D_1*w25 + D_2*w24 + tmp10;
2148                              const double tmp15_1 = D_2*w55;                              EM_S[INDEX2(2,3,4)]+=tmp3 + tmp4;
2149                              const double tmp16_1 = D_1*w56;                              EM_S[INDEX2(3,0,4)]+=tmp2;
2150                              const double tmp17_1 = D_1*w55;                              EM_S[INDEX2(3,1,4)]+=tmp6 + tmp7;
2151                              const double tmp18_1 = D_2*w56;                              EM_S[INDEX2(3,2,4)]+=tmp3 + tmp4;
2152                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(3,3,4)]+=D_0*w25 + D_3*w24 + tmp5;
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp2_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp2_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp2_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;  
2153                          } else { // constant data                          } else { // constant data
2154                              const double tmp0_1 = D_p[0]*w57;                              const double D_0 = D_p[0];
2155                              const double tmp1_1 = D_p[0]*w58;                              EM_S[INDEX2(0,0,4)]+=16*D_0*w23;
2156                              const double tmp2_1 = D_p[0]*w59;                              EM_S[INDEX2(0,1,4)]+=8*D_0*w23;
2157                              EM_S[INDEX2(0,0,4)]+=tmp2_1;                              EM_S[INDEX2(0,2,4)]+=8*D_0*w23;
2158                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(0,3,4)]+=4*D_0*w23;
2159                              EM_S[INDEX2(2,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=8*D_0*w23;
2160                              EM_S[INDEX2(3,0,4)]+=tmp1_1;                              EM_S[INDEX2(1,1,4)]+=16*D_0*w23;
2161                              EM_S[INDEX2(0,1,4)]+=tmp0_1;                              EM_S[INDEX2(1,2,4)]+=4*D_0*w23;
2162                              EM_S[INDEX2(1,1,4)]+=tmp2_1;                              EM_S[INDEX2(1,3,4)]+=8*D_0*w23;
2163                              EM_S[INDEX2(2,1,4)]+=tmp1_1;                              EM_S[INDEX2(2,0,4)]+=8*D_0*w23;
2164                              EM_S[INDEX2(3,1,4)]+=tmp0_1;                              EM_S[INDEX2(2,1,4)]+=4*D_0*w23;
2165                              EM_S[INDEX2(0,2,4)]+=tmp0_1;                              EM_S[INDEX2(2,2,4)]+=16*D_0*w23;
2166                              EM_S[INDEX2(1,2,4)]+=tmp1_1;                              EM_S[INDEX2(2,3,4)]+=8*D_0*w23;
2167                              EM_S[INDEX2(2,2,4)]+=tmp2_1;                              EM_S[INDEX2(3,0,4)]+=4*D_0*w23;
2168                              EM_S[INDEX2(3,2,4)]+=tmp0_1;                              EM_S[INDEX2(3,1,4)]+=8*D_0*w23;
2169                              EM_S[INDEX2(0,3,4)]+=tmp1_1;                              EM_S[INDEX2(3,2,4)]+=8*D_0*w23;
2170                              EM_S[INDEX2(1,3,4)]+=tmp0_1;                              EM_S[INDEX2(3,3,4)]+=16*D_0*w23;
                             EM_S[INDEX2(2,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp2_1;  
2171                          }                          }
2172                      }                      }
2173                      ///////////////                      ///////////////
# Line 2090  void Rectangle::assemblePDESingle(Paso_S Line 2185  void Rectangle::assemblePDESingle(Paso_S
2185                              const double X_1_2 = X_p[INDEX2(1,2,2)];                              const double X_1_2 = X_p[INDEX2(1,2,2)];
2186                              const double X_0_3 = X_p[INDEX2(0,3,2)];                              const double X_0_3 = X_p[INDEX2(0,3,2)];
2187                              const double X_1_3 = X_p[INDEX2(1,3,2)];                              const double X_1_3 = X_p[INDEX2(1,3,2)];
2188                              const double tmp0_0 = X_0_2 + X_0_3;                              const double tmp0 = 6*w15*(X_1_1 + X_1_3);
2189                              const double tmp1_0 = X_1_1 + X_1_3;                              const double tmp1 = 6*w16*(X_0_2 + X_0_3);
2190                              const double tmp2_0 = X_1_0 + X_1_2;                              const double tmp2 = 6*w11*(X_0_0 + X_0_1);
2191                              const double tmp3_0 = X_0_0 + X_0_1;                              const double tmp3 = 6*w12*(X_1_0 + X_1_2);
2192                              const double tmp0_1 = tmp0_0*w63;                              const double tmp4 = 6*w15*(X_1_0 + X_1_2);
2193                              const double tmp1_1 = tmp1_0*w62;                              const double tmp5 = w27*(X_0_2 + X_0_3);
2194                              const double tmp2_1 = tmp2_0*w61;                              const double tmp6 = w26*(X_0_0 + X_0_1);
2195                              const double tmp3_1 = tmp3_0*w60;                              const double tmp7 = 6*w12*(X_1_1 + X_1_3);
2196                              const double tmp4_1 = tmp0_0*w65;                              const double tmp8 = w29*(X_1_1 + X_1_3);
2197                              const double tmp5_1 = tmp3_0*w64;                              const double tmp9 = w27*(-X_0_0 - X_0_1);
2198                              const double tmp6_1 = tmp2_0*w62;                              const double tmp10 = w28*(X_1_0 + X_1_2);
2199                              const double tmp7_1 = tmp1_0*w61;                              const double tmp11 = w26*(-X_0_2 - X_0_3);
2200                              const double tmp8_1 = tmp2_0*w66;                              const double tmp12 = w29*(X_1_0 + X_1_2);
2201                              const double tmp9_1 = tmp3_0*w63;                              const double tmp13 = w27*(X_0_0 + X_0_1);
2202                              const double tmp10_1 = tmp0_0*w60;                              const double tmp14 = w28*(X_1_1 + X_1_3);
2203                              const double tmp11_1 = tmp1_0*w67;                              const double tmp15 = w26*(X_0_2 + X_0_3);
2204                              const double tmp12_1 = tmp1_0*w66;                              EM_F[0]+=tmp0 + tmp1 + tmp2 + tmp3;
2205                              const double tmp13_1 = tmp3_0*w65;                              EM_F[1]+=tmp4 + tmp5 + tmp6 + tmp7;
2206                              const double tmp14_1 = tmp0_0*w64;                              EM_F[2]+=tmp10 + tmp11 + tmp8 + tmp9;
2207                              const double tmp15_1 = tmp2_0*w67;                              EM_F[3]+=tmp12 + tmp13 + tmp14 + tmp15;
                             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;  
2208                          } else { // constant data                          } else { // constant data
2209                              const double tmp0_1 = X_p[1]*w69;                              const double X_0 = X_p[0];
2210                              const double tmp1_1 = X_p[0]*w68;                              const double X_1 = X_p[1];
2211                              const double tmp2_1 = X_p[0]*w70;                              EM_F[0]+=3*X_0*w19 + 6*X_1*w20;
2212                              const double tmp3_1 = X_p[1]*w71;                              EM_F[1]+=-3*X_0*w19 + 6*X_1*w20;
2213                              EM_F[0]+=tmp0_1 + tmp1_1;                              EM_F[2]+=3*X_0*w19 - 6*X_1*w20;
2214                              EM_F[1]+=tmp0_1 + tmp2_1;                              EM_F[3]+=-3*X_0*w19 - 6*X_1*w20;
                             EM_F[2]+=tmp1_1 + tmp3_1;  
                             EM_F[3]+=tmp2_1 + tmp3_1;  
2215                          }                          }
2216                      }                      }
2217                      ///////////////                      ///////////////
# Line 2136  void Rectangle::assemblePDESingle(Paso_S Line 2225  void Rectangle::assemblePDESingle(Paso_S
2225                              const double Y_1 = Y_p[1];                              const double Y_1 = Y_p[1];
2226                              const double Y_2 = Y_p[2];                              const double Y_2 = Y_p[2];
2227                              const double Y_3 = Y_p[3];                              const double Y_3 = Y_p[3];
2228                              const double tmp0_0 = Y_1 + Y_2;                              const double tmp0 = 6*w23*(Y_1 + Y_2);
2229                              const double tmp1_0 = Y_0 + Y_3;                              const double tmp1 = 6*w23*(Y_0 + Y_3);
2230                              const double tmp0_1 = Y_0*w72;                              EM_F[0]+=6*Y_0*w21 + 6*Y_3*w22 + tmp0;
2231                              const double tmp1_1 = tmp0_0*w73;                              EM_F[1]+=6*Y_1*w21 + 6*Y_2*w22 + tmp1;
2232                              const double tmp2_1 = Y_3*w74;                              EM_F[2]+=6*Y_1*w22 + 6*Y_2*w21 + tmp1;
2233                              const double tmp3_1 = Y_1*w72;                              EM_F[3]+=6*Y_0*w22 + 6*Y_3*w21 + tmp0;
                             const double tmp4_1 = tmp1_0*w73;  
                             const double tmp5_1 = Y_2*w74;  
                             const double tmp6_1 = Y_2*w72;  
                             const double tmp7_1 = Y_1*w74;  
                             const double tmp8_1 = Y_3*w72;  
                             const double tmp9_1 = Y_0*w74;  
                             EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;  
                             EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;  
                             EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;  
2234                          } else { // constant data                          } else { // constant data
2235                              const double tmp0_1 = Y_p[0]*w75;                              const double Y_0 = Y_p[0];
2236                              EM_F[0]+=tmp0_1;                              EM_F[0]+=36*Y_0*w23;
2237                              EM_F[1]+=tmp0_1;                              EM_F[1]+=36*Y_0*w23;
2238                              EM_F[2]+=tmp0_1;                              EM_F[2]+=36*Y_0*w23;
2239                              EM_F[3]+=tmp0_1;                              EM_F[3]+=36*Y_0*w23;
2240                          }                          }
2241                      }                      }
2242                        /* GENERATOR SNIP_PDE_SINGLE BOTTOM */
2243    
2244                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2245                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
2246                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2247                  } // end k0 loop                  } // end k0 loop
2248              } // end k1 loop              } // end k1 loop
# Line 2176  void Rectangle::assemblePDESingleReduced Line 2256  void Rectangle::assemblePDESingleReduced
2256          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2257          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
2258  {  {
2259      const double h0 = m_l0/m_gNE0;      const double w0 = .25;
2260      const double h1 = m_l1/m_gNE1;      const double w1 = .125*m_dx[0];
2261      const double w0 = -.25*h1/h0;      const double w2 = .125*m_dx[1];
2262      const double w1 = .25;      const double w3 = .0625*m_dx[0]*m_dx[1];
2263      const double w2 = -.25;      const double w4 = .25*m_dx[0]/m_dx[1];
2264      const double w3 = .25*h0/h1;      const double w5 = .25*m_dx[1]/m_dx[0];
     const double w4 = -.25*h0/h1;  
     const double w5 = .25*h1/h0;  
     const double w6 = -.125*h1;  
     const double w7 = -.125*h0;  
     const double w8 = .125*h1;  
     const double w9 = .125*h0;  
     const double w10 = .0625*h0*h1;  
     const double w11 = -.5*h1;  
     const double w12 = -.5*h0;  
     const double w13 = .5*h1;  
     const double w14 = .5*h0;  
     const double w15 = .25*h0*h1;  
2265    
2266      rhs.requireWrite();      rhs.requireWrite();
2267  #pragma omp parallel  #pragma omp parallel
2268      {      {
2269          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2270  #pragma omp for  #pragma omp for
2271              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {              for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
2272                  for (index_t k0=0; k0<m_NE0; ++k0)  {                  for (index_t k0=0; k0<m_NE[0]; ++k0)  {
2273                      bool add_EM_S=false;                      bool add_EM_S=false;
2274                      bool add_EM_F=false;                      bool add_EM_F=false;
2275                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
2276                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
2277                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE[0]*k1;
2278                      ///////////////                      ///////////////
2279                      // process A //                      // process A //
2280                      ///////////////                      ///////////////
# Line 2217  void Rectangle::assemblePDESingleReduced Line 2285  void Rectangle::assemblePDESingleReduced
2285                          const double A_10 = A_p[INDEX2(1,0,2)];                          const double A_10 = A_p[INDEX2(1,0,2)];
2286                          const double A_01 = A_p[INDEX2(0,1,2)];                          const double A_01 = A_p[INDEX2(0,1,2)];
2287                          const double A_11 = A_p[INDEX2(1,1,2)];                          const double A_11 = A_p[INDEX2(1,1,2)];
2288                          const double tmp0_0 = A_01 + A_10;                          const double tmp0 = (A_01 + A_10)*w0;
2289                          const double tmp0_1 = A_11*w3;                          const double tmp1 = A_00*w5;
2290                          const double tmp1_1 = A_00*w0;                          const double tmp2 = A_01*w0;
2291                          const double tmp2_1 = A_01*w1;                          const double tmp3 = A_10*w0;
2292                          const double tmp3_1 = A_10*w2;                          const double tmp4 = A_11*w4;
2293                          const double tmp4_1 = tmp0_0*w1;                          EM_S[INDEX2(0,0,4)]+=tmp4 + tmp0 + tmp1;
2294                          const double tmp5_1 = A_11*w4;                          EM_S[INDEX2(1,0,4)]+=tmp4 - tmp1 + tmp3 - tmp2;
2295                          const double tmp6_1 = A_00*w5;                          EM_S[INDEX2(2,0,4)]+=tmp2 - tmp3 - tmp4 + tmp1;
2296                          const double tmp7_1 = tmp0_0*w2;                          EM_S[INDEX2(3,0,4)]+=-tmp1 - tmp4 - tmp0;
2297                          const double tmp8_1 = A_10*w1;                          EM_S[INDEX2(0,1,4)]+=tmp4 - tmp1 + tmp2 - tmp3;
2298                          const double tmp9_1 = A_01*w2;                          EM_S[INDEX2(1,1,4)]+=tmp4 + tmp1 - tmp0;
2299                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                          EM_S[INDEX2(2,1,4)]+=-tmp1 + tmp0 - tmp4;
2300                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(3,1,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;
2301                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                          EM_S[INDEX2(0,2,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;
2302                          EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;                          EM_S[INDEX2(1,2,4)]+=-tmp1 + tmp0 - tmp4;
2303                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                          EM_S[INDEX2(2,2,4)]+=tmp4 + tmp1 - tmp0;
2304                          EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;                          EM_S[INDEX2(3,2,4)]+=tmp4 - tmp1 + tmp2 - tmp3;
2305                          EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;                          EM_S[INDEX2(0,3,4)]+=-tmp1 - tmp4 - tmp0;
2306                          EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(1,3,4)]+=tmp2 - tmp3 - tmp4 + tmp1;
2307                          EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(2,3,4)]+=tmp4 - tmp1 + tmp3 - tmp2;
2308                          EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;                          EM_S[INDEX2(3,3,4)]+=tmp4 + tmp0 + tmp1;
                         EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;  
2309                      }                      }
2310                      ///////////////                      ///////////////
2311                      // process B //                      // process B //
# Line 2251  void Rectangle::assemblePDESingleReduced Line 2313  void Rectangle::assemblePDESingleReduced
2313                      if (!B.isEmpty()) {                      if (!B.isEmpty()) {
2314                          add_EM_S=true;                          add_EM_S=true;
2315                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2316                          const double tmp2_1 = B_p[0]*w8;                          const double tmp0 = B_p[0]*w2;
2317                          const double tmp0_1 = B_p[0]*w6;                          const double tmp1 = B_p[1]*w1;
2318                          const double tmp3_1 = B_p[1]*w9;                          EM_S[INDEX2(0,0,4)]+=-tmp0 - tmp1;
2319                          const double tmp1_1 = B_p[1]*w7;                          EM_S[INDEX2(1,0,4)]+= tmp0 - tmp1;
2320                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,0,4)]+= tmp1 - tmp0;
2321                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(3,0,4)]+= tmp0 + tmp1;
2322                          EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(0,1,4)]+=-tmp0 - tmp1;
2323                          EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(1,1,4)]+= tmp0 - tmp1;
2324                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,1,4)]+= tmp1 - tmp0;
2325                          EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(3,1,4)]+= tmp0 + tmp1;
2326                          EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(0,2,4)]+=-tmp0 - tmp1;
2327                          EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(1,2,4)]+= tmp0 - tmp1;
2328                          EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,2,4)]+= tmp1 - tmp0;
2329                          EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(3,2,4)]+= tmp0 + tmp1;
2330                          EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(0,3,4)]+=-tmp0 - tmp1;
2331                          EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(1,3,4)]+= tmp0 - tmp1;
2332                          EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,3,4)]+= tmp1 - tmp0;
2333                          EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(3,3,4)]+= tmp0 + tmp1;
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;  
2334                      }                      }
2335                      ///////////////                      ///////////////
2336                      // process C //                      // process C //
# Line 2278  void Rectangle::assemblePDESingleReduced Line 2338  void Rectangle::assemblePDESingleReduced
2338                      if (!C.isEmpty()) {                      if (!C.isEmpty()) {
2339                          add_EM_S=true;                          add_EM_S=true;
2340                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2341                          const double tmp1_1 = C_p[1]*w7;                          const double tmp0 = C_p[0]*w2;
2342                          const double tmp0_1 = C_p[0]*w8;                          const double tmp1 = C_p[1]*w1;
2343                          const double tmp3_1 = C_p[0]*w6;                          EM_S[INDEX2(0,0,4)]+=-tmp1 - tmp0;
2344                          const double tmp2_1 = C_p[1]*w9;                          EM_S[INDEX2(1,0,4)]+=-tmp1 - tmp0;
2345                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(2,0,4)]+=-tmp1 - tmp0;
2346                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(3,0,4)]+=-tmp1 - tmp0;
2347                          EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(0,1,4)]+= tmp0 - tmp1;
2348                          EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(1,1,4)]+= tmp0 - tmp1;
2349                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,1,4)]+= tmp0 - tmp1;
2350                          EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(3,1,4)]+= tmp0 - tmp1;
2351                          EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(0,2,4)]+= tmp1 - tmp0;
2352                          EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(1,2,4)]+= tmp1 - tmp0;
2353                          EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(2,2,4)]+= tmp1 - tmp0;
2354                          EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(3,2,4)]+= tmp1 - tmp0;
2355                          EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(0,3,4)]+= tmp0 + tmp1;
2356                          EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(1,3,4)]+= tmp0 + tmp1;
2357                          EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;                          EM_S[INDEX2(2,3,4)]+= tmp0 + tmp1;
2358                          EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;                          EM_S[INDEX2(3,3,4)]+= tmp0 + tmp1;
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;  
2359                      }                      }
2360                      ///////////////                      ///////////////
2361                      // process D //                      // process D //
# Line 2305  void Rectangle::assemblePDESingleReduced Line 2363  void Rectangle::assemblePDESingleReduced
2363                      if (!D.isEmpty()) {                      if (!D.isEmpty()) {
2364                          add_EM_S=true;                          add_EM_S=true;
2365                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2366                          const double tmp0_1 = D_p[0]*w10;                          EM_S[INDEX2(0,0,4)]+=D_p[0]*w3;
2367                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(1,0,4)]+=D_p[0]*w3;
2368                          EM_S[INDEX2(1,0,4)]+=tmp0_1;                          EM_S[INDEX2(2,0,4)]+=D_p[0]*w3;
2369                          EM_S[INDEX2(2,0,4)]+=tmp0_1;                          EM_S[INDEX2(3,0,4)]+=D_p[0]*w3;
2370                          EM_S[INDEX2(3,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,1,4)]+=D_p[0]*w3;
2371                          EM_S[INDEX2(0,1,4)]+=tmp0_1;                          EM_S[INDEX2(1,1,4)]+=D_p[0]*w3;
2372                          EM_S[INDEX2(1,1,4)]+=tmp0_1;                          EM_S[INDEX2(2,1,4)]+=D_p[0]*w3;
2373                          EM_S[INDEX2(2,1,4)]+=tmp0_1;                          EM_S[INDEX2(3,1,4)]+=D_p[0]*w3;
2374                          EM_S[INDEX2(3,1,4)]+=tmp0_1;                          EM_S[INDEX2(0,2,4)]+=D_p[0]*w3;
2375                          EM_S[INDEX2(0,2,4)]+=tmp0_1;                          EM_S[INDEX2(1,2,4)]+=D_p[0]*w3;
2376                          EM_S[INDEX2(1,2,4)]+=tmp0_1;                          EM_S[INDEX2(2,2,4)]+=D_p[0]*w3;
2377                          EM_S[INDEX2(2,2,4)]+=tmp0_1;                          EM_S[INDEX2(3,2,4)]+=D_p[0]*w3;
2378                          EM_S[INDEX2(3,2,4)]+=tmp0_1;                          EM_S[INDEX2(0,3,4)]+=D_p[0]*w3;
2379                          EM_S[INDEX2(0,3,4)]+=tmp0_1;                          EM_S[INDEX2(1,3,4)]+=D_p[0]*w3;
2380                          EM_S[INDEX2(1,3,4)]+=tmp0_1;                          EM_S[INDEX2(2,3,4)]+=D_p[0]*w3;
2381                          EM_S[INDEX2(2,3,4)]+=tmp0_1;                          EM_S[INDEX2(3,3,4)]+=D_p[0]*w3;
                         EM_S[INDEX2(3,3,4)]+=tmp0_1;  
2382                      }                      }
2383                      ///////////////                      ///////////////
2384                      // process X //                      // process X //
# Line 2329  void Rectangle::assemblePDESingleReduced Line 2386  void Rectangle::assemblePDESingleReduced
2386                      if (!X.isEmpty()) {                      if (!X.isEmpty()) {
2387                          add_EM_F=true;                          add_EM_F=true;
2388                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2389                          const double tmp0_1 = X_p[0]*w11;                          const double tmp0 = 4*X_p[0]*w2;
2390                          const double tmp2_1 = X_p[0]*w13;                          const double tmp1 = 4*X_p[1]*w1;
2391                          const double tmp1_1 = X_p[1]*w12;                          EM_F[0]+=-tmp0 - tmp1;
2392                          const double tmp3_1 = X_p[1]*w14;                          EM_F[1]+=-tmp1 + tmp0;
2393                          EM_F[0]+=tmp0_1 + tmp1_1;                          EM_F[2]+=-tmp0 + tmp1;
2394                          EM_F[1]+=tmp1_1 + tmp2_1;                          EM_F[3]+= tmp0 + tmp1;
                         EM_F[2]+=tmp0_1 + tmp3_1;  
                         EM_F[3]+=tmp2_1 + tmp3_1;  
2395                      }                      }
2396                      ///////////////                      ///////////////
2397                      // process Y //                      // process Y //
# Line 2344  void Rectangle::assemblePDESingleReduced Line 2399  void Rectangle::assemblePDESingleReduced
2399                      if (!Y.isEmpty()) {                      if (!Y.isEmpty()) {
2400                          add_EM_F=true;                          add_EM_F=true;
2401                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2402                          const double tmp0_1 = Y_p[0]*w15;                          EM_F[0]+=4*Y_p[0]*w3;
2403                          EM_F[0]+=tmp0_1;                          EM_F[1]+=4*Y_p[0]*w3;
2404                          EM_F[1]+=tmp0_1;                          EM_F[2]+=4*Y_p[0]*w3;
2405                          EM_F[2]+=tmp0_1;                          EM_F[3]+=4*Y_p[0]*w3;
                         EM_F[3]+=tmp0_1;  
2406                      }                      }
2407    
2408                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2409                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
2410                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2411                  } // end k0 loop                  } // end k0 loop
2412              } // end k1 loop              } // end k1 loop
# Line 2366  void Rectangle::assemblePDESystem(Paso_S Line 2420  void Rectangle::assemblePDESystem(Paso_S
2420          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2421          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
2422  {  {
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
2423      dim_t numEq, numComp;      dim_t numEq, numComp;
2424      if (!mat)      if (!mat)
2425          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
# Line 2376  void Rectangle::assemblePDESystem(Paso_S Line 2428  void Rectangle::assemblePDESystem(Paso_S
2428          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
2429      }      }
2430    
2431      const double w0 = -0.1555021169820365539*h1/h0;      /* GENERATOR SNIP_PDE_SYSTEM_PRE TOP */
2432        const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];
2433      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
2434      const double w2 = -0.15550211698203655390;      const double w2 = -0.15550211698203655390;
2435      const double w3 = 0.041666666666666666667*h0/h1;      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];
2436      const double w4 = 0.15550211698203655390;      const double w4 = -0.01116454968463011277*m_dx[1]/m_dx[0];
2437      const double w5 = -0.041666666666666666667;      const double w5 = 0.011164549684630112770;
2438      const double w6 = -0.01116454968463011277*h1/h0;      const double w6 = -0.041666666666666666667*m_dx[1]/m_dx[0];
2439      const double w7 = 0.011164549684630112770;      const double w7 = 0.1555021169820365539*m_dx[0]/m_dx[1];
2440      const double w8 = -0.011164549684630112770;      const double w8 = 0.01116454968463011277*m_dx[0]/m_dx[1];
2441      const double w9 = -0.041666666666666666667*h1/h0;      const double w9 = -0.25;
2442      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.16666666666666666667*m_dx[0]/m_dx[1];
2443      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = -0.032861463941450536761*m_dx[1];
2444      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = -0.032861463941450536761*m_dx[0];
2445      const double w13 = 0.01116454968463011277*h0/h1;      const double w13 = -0.12264065304058601714*m_dx[1];
2446      const double w14 = 0.01116454968463011277*h1/h0;      const double w14 = -0.0023593469594139828636*m_dx[1];
2447      const double w15 = 0.041666666666666666667*h1/h0;      const double w15 = -0.008805202725216129906*m_dx[0];
2448      const double w16 = -0.01116454968463011277*h0/h1;      const double w16 = -0.008805202725216129906*m_dx[1];
2449      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.12264065304058601714*m_dx[0];
2450      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.0023593469594139828636*m_dx[0];
2451      const double w19 = 0.25000000000000000000;      const double w19 = -0.16666666666666666667*m_dx[1];
2452      const double w20 = -0.25000000000000000000;      const double w20 = -0.083333333333333333333*m_dx[0];
2453      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.025917019497006092316*m_dx[0]*m_dx[1];
2454      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = 0.0018607582807716854616*m_dx[0]*m_dx[1];
2455      const double w23 = -0.16666666666666666667*h0/h1;      const double w23 = 0.0069444444444444444444*m_dx[0]*m_dx[1];
2456      const double w24 = 0.33333333333333333333*h1/h0;      const double w24 = 0.09672363354357992482*m_dx[0]*m_dx[1];
2457      const double w25 = 0.33333333333333333333*h0/h1;      const double w25 = 0.000498588678642297402*m_dx[0]*m_dx[1];
2458      const double w26 = 0.16666666666666666667*h1/h0;      const double w26 = 0.19716878364870322056*m_dx[1];
2459      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = 0.052831216351296779436*m_dx[1];
2460      const double w28 = -0.032861463941450536761*h1;      const double w28 = 0.19716878364870322056*m_dx[0];
2461      const double w29 = -0.032861463941450536761*h0;      const double w29 = 0.052831216351296779436*m_dx[0];
2462      const double w30 = -0.12264065304058601714*h1;      /* GENERATOR SNIP_PDE_SYSTEM_PRE BOTTOM */
     const double w31 = -0.0023593469594139828636*h1;  
     const double w32 = -0.008805202725216129906*h0;  
     const double w33 = -0.008805202725216129906*h1;  
     const double w34 = 0.032861463941450536761*h1;  
     const double w35 = 0.008805202725216129906*h1;  
     const double w36 = 0.008805202725216129906*h0;  
     const double w37 = 0.0023593469594139828636*h1;  
     const double w38 = 0.12264065304058601714*h1;  
     const double w39 = 0.032861463941450536761*h0;  
     const double w40 = -0.12264065304058601714*h0;  
     const double w41 = -0.0023593469594139828636*h0;  
     const double w42 = 0.0023593469594139828636*h0;  
     const double w43 = 0.12264065304058601714*h0;  
     const double w44 = -0.16666666666666666667*h1;  
     const double w45 = -0.083333333333333333333*h0;  
     const double w46 = 0.083333333333333333333*h1;  
     const double w47 = 0.16666666666666666667*h1;  
     const double w48 = 0.083333333333333333333*h0;  
     const double w49 = -0.16666666666666666667*h0;  
     const double w50 = 0.16666666666666666667*h0;  
     const double w51 = -0.083333333333333333333*h1;  
     const double w52 = 0.025917019497006092316*h0*h1;  
     const double w53 = 0.0018607582807716854616*h0*h1;  
     const double w54 = 0.0069444444444444444444*h0*h1;  
     const double w55 = 0.09672363354357992482*h0*h1;  
     const double w56 = 0.00049858867864229740201*h0*h1;  
     const double w57 = 0.055555555555555555556*h0*h1;  
     const double w58 = 0.027777777777777777778*h0*h1;  
     const double w59 = 0.11111111111111111111*h0*h1;  
     const double w60 = -0.19716878364870322056*h1;  
     const double w61 = -0.19716878364870322056*h0;  
     const double w62 = -0.052831216351296779436*h0;  
     const double w63 = -0.052831216351296779436*h1;  
     const double w64 = 0.19716878364870322056*h1;  
     const double w65 = 0.052831216351296779436*h1;  
     const double w66 = 0.19716878364870322056*h0;  
     const double w67 = 0.052831216351296779436*h0;  
     const double w68 = -0.5*h1;  
     const double w69 = -0.5*h0;  
     const double w70 = 0.5*h1;  
     const double w71 = 0.5*h0;  
     const double w72 = 0.1555021169820365539*h0*h1;  
     const double w73 = 0.041666666666666666667*h0*h1;  
     const double w74 = 0.01116454968463011277*h0*h1;  
     const double w75 = 0.25*h0*h1;  
2463    
2464      rhs.requireWrite();      rhs.requireWrite();
2465  #pragma omp parallel  #pragma omp parallel
2466      {      {
2467          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2468  #pragma omp for  #pragma omp for
2469              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {              for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
2470                  for (index_t k0=0; k0<m_NE0; ++k0)  {                  for (index_t k0=0; k0<m_NE[0]; ++k0)  {
2471                      bool add_EM_S=false;                      bool add_EM_S=false;
2472                      bool add_EM_F=false;                      bool add_EM_F=false;
2473                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
2474                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
2475                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE[0]*k1;
2476                        /* GENERATOR SNIP_PDE_SYSTEM TOP */
2477                      ///////////////                      ///////////////
2478                      // process A //                      // process A //
2479                      ///////////////                      ///////////////
2480                      if (!A.isEmpty()) {                      if (!A.isEmpty()) {
2481                          add_EM_S=true;                          add_EM_S = true;
2482                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p = const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2483                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
2484                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2485                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2486                                      const double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];                                      const double A_00_0 = A_p[INDEX5(k,0,m,0,0,numEq,2,numComp,2)];
2487                                      const double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];                                      const double A_01_0 = A_p[INDEX5(k,0,m,1,0,numEq,2,numComp,2)];
2488                                      const double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];                                      const double A_10_0 = A_p[INDEX5(k,1,m,0,0,numEq,2,numComp,2)];
2489                                      const double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];                                      const double A_11_0 = A_p[INDEX5(k,1,m,1,0,numEq,2,numComp,2)];
2490                                      const double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];                                      const double A_00_1 = A_p[INDEX5(k,0,m,0,1,numEq,2,numComp,2)];
2491                                      const double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];                                      const double A_01_1 = A_p[INDEX5(k,0,m,1,1,numEq,2,numComp,2)];
2492                                      const double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];                                      const double A_10_1 = A_p[INDEX5(k,1,m,0,1,numEq,2,numComp,2)];
2493                                      const double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];                                      const double A_11_1 = A_p[INDEX5(k,1,m,1,1,numEq,2,numComp,2)];
2494                                      const double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];                                      const double A_00_2 = A_p[INDEX5(k,0,m,0,2,numEq,2,numComp,2)];
2495                                      const double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];                                      const double A_01_2 = A_p[INDEX5(k,0,m,1,2,numEq,2,numComp,2)];
2496                                      const double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];                                      const double A_10_2 = A_p[INDEX5(k,1,m,0,2,numEq,2,numComp,2)];
2497                                      const double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];                                      const double A_11_2 = A_p[INDEX5(k,1,m,1,2,numEq,2,numComp,2)];
2498                                      const double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];                                      const double A_00_3 = A_p[INDEX5(k,0,m,0,3,numEq,2,numComp,2)];
2499                                      const double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];                                      const double A_01_3 = A_p[INDEX5(k,0,m,1,3,numEq,2,numComp,2)];
2500                                      const double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];                                      const double A_10_3 = A_p[INDEX5(k,1,m,0,3,numEq,2,numComp,2)];
2501                                      const double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];                                      const double A_11_3 = A_p[INDEX5(k,1,m,1,3,numEq,2,numComp,2)];
2502                                      const double tmp0_0 = A_01_0 + A_01_3;                                      const double tmp0 = w3*(A_11_0 + A_11_1 + A_11_2 + A_11_3);
2503                                      const double tmp1_0 = A_00_0 + A_00_1;                                      const double tmp1 = w1*(A_01_0 + A_01_3 - A_10_1 - A_10_2);
2504                                      const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                                      const double tmp2 = w4*(A_00_2 + A_00_3);
2505                                      const double tmp3_0 = A_00_2 + A_00_3;                                      const double tmp3 = w0*(A_00_0 + A_00_1);
2506                                      const double tmp4_0 = A_10_1 + A_10_2;                                      const double tmp4 = w5*(A_01_2 - A_10_3);
2507                                      const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                                      const double tmp5 = w2*(-A_01_1 + A_10_0);
2508                                      const double tmp6_0 = A_01_3 + A_10_0;                                      const double tmp6 = w5*(A_01_3 + A_10_0);
2509                                      const double tmp7_0 = A_01_0 + A_10_3;                                      const double tmp7 = w3*(-A_11_0 - A_11_1 - A_11_2 - A_11_3);
2510                                      const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                                      const double tmp8 = w6*(A_00_0 + A_00_1 + A_00_2 + A_00_3);
2511                                      const double tmp9_0 = A_01_0 + A_10_0;                                      const double tmp9 = w1*(A_01_1 + A_01_2 + A_10_1 + A_10_2);
2512                                      const double tmp10_0 = A_01_3 + A_10_3;                                      const double tmp10 = w2*(-A_01_0 - A_10_3);
2513                                      const double tmp11_0 = A_11_1 + A_11_3;                                      const double tmp11 = w4*(A_00_0 + A_00_1);
2514                                      const double tmp12_0 = A_11_0 + A_11_2;                                      const double tmp12 = w0*(A_00_2 + A_00_3);
2515                                      const double tmp13_0 = A_01_2 + A_10_1;                                      const double tmp13 = w5*(A_01_1 - A_10_0);
2516                                      const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                                      const double tmp14 = w2*(-A_01_2 + A_10_3);
2517                                      const double tmp15_0 = A_01_1 + A_10_2;                                      const double tmp15 = w7*(A_11_0 + A_11_2);
2518                                      const double tmp16_0 = A_10_0 + A_10_3;                                      const double tmp16 = w4*(-A_00_2 - A_00_3);
2519                                      const double tmp17_0 = A_01_1 + A_01_2;                                      const double tmp17 = w0*(-A_00_0 - A_00_1);
2520                                      const double tmp18_0 = A_01_1 + A_10_1;                                      const double tmp18 = w5*(A_01_3 + A_10_3);
2521                                      const double tmp19_0 = A_01_2 + A_10_2;                                      const double tmp19 = w8*(A_11_1 + A_11_3);
2522                                      const double tmp0_1 = A_10_3*w8;                                      const double tmp20 = w2*(-A_01_0 - A_10_0);
2523                                      const double tmp1_1 = tmp0_0*w1;                                      const double tmp21 = w7*(A_11_1 + A_11_3);
2524                                      const double tmp2_1 = A_01_1*w4;                                      const double tmp22 = w4*(-A_00_0 - A_00_1);
2525                                      const double tmp3_1 = tmp1_0*w0;                                      const double tmp23 = w0*(-A_00_2 - A_00_3);
2526                                      const double tmp4_1 = A_01_2*w7;                                      const double tmp24 = w5*(A_01_0 + A_10_0);
2527                                      const double tmp5_1 = tmp2_0*w3;                                      const double tmp25 = w8*(A_11_0 + A_11_2);
2528                                      const double tmp6_1 = tmp3_0*w6;                                      const double tmp26 = w2*(-A_01_3 - A_10_3);
2529                                      const double tmp7_1 = A_10_0*w2;                                      const double tmp27 = w5*(-A_01_1 - A_10_2);
2530                                      const double tmp8_1 = tmp4_0*w5;                                      const double tmp28 = w1*(-A_01_0 - A_01_3 - A_10_0 - A_10_3);
2531                                      const double tmp9_1 = tmp2_0*w10;                                      const double tmp29 = w2*(A_01_2 + A_10_1);
2532                                      const double tmp10_1 = tmp5_0*w9;                                      const double tmp30 = w7*(-A_11_1 - A_11_3);
2533                                      const double tmp11_1 = tmp6_0*w7;                                      const double tmp31 = w1*(-A_01_1 - A_01_2 + A_10_0 + A_10_3);
2534                                      const double tmp12_1 = tmp7_0*w4;                                      const double tmp32 = w5*(-A_01_0 + A_10_2);
2535                                      const double tmp13_1 = tmp8_0*w1;                                      const double tmp33 = w8*(-A_11_0 - A_11_2);
2536                                      const double tmp14_1 = A_10_0*w8;                                      const double tmp34 = w6*(-A_00_0 - A_00_1 - A_00_2 - A_00_3);
2537                                      const double tmp15_1 = A_01_2*w4;                                      const double tmp35 = w2*(A_01_3 - A_10_1);
2538                                      const double tmp16_1 = tmp3_0*w0;                                      const double tmp36 = w5*(A_01_0 + A_10_3);
2539                                      const double tmp17_1 = A_01_1*w7;                                      const double tmp37 = w2*(-A_01_3 - A_10_0);
2540                                      const double tmp18_1 = tmp1_0*w6;                                      const double tmp38 = w7*(-A_11_0 - A_11_2);
2541                                      const double tmp19_1 = A_10_3*w2;                                      const double tmp39 = w5*(-A_01_3 + A_10_1);
2542                                      const double tmp20_1 = tmp9_0*w4;                                      const double tmp40 = w8*(-A_11_1 - A_11_3);
2543                                      const double tmp21_1 = tmp1_0*w11;                                      const double tmp41 = w2*(A_01_0 - A_10_2);
2544                                      const double tmp22_1 = tmp10_0*w7;                                      const double tmp42 = w5*(A_01_1 - A_10_3);
2545                                      const double tmp23_1 = tmp3_0*w14;                                      const double tmp43 = w2*(-A_01_2 + A_10_0);
2546                                      const double tmp24_1 = tmp11_0*w13;                                      const double tmp44 = w5*(A_01_2 - A_10_0);
2547                                      const double tmp25_1 = tmp12_0*w12;                                      const double tmp45 = w2*(-A_01_1 + A_10_3);
2548                                      const double tmp26_1 = tmp10_0*w4;                                      const double tmp46 = w5*(-A_01_0 + A_10_1);
2549                                      const double tmp27_1 = tmp3_0*w11;                                      const double tmp47 = w2*(A_01_3 - A_10_2);
2550                                      const double tmp28_1 = tmp9_0*w7;                                      const double tmp48 = w5*(-A_01_1 - A_10_1);
2551                                      const double tmp29_1 = tmp1_0*w14;                                      const double tmp49 = w2*(A_01_2 + A_10_2);
2552                                      const double tmp30_1 = tmp12_0*w13;                                      const double tmp50 = w5*(-A_01_3 + A_10_2);
2553                                      const double tmp31_1 = tmp11_0*w12;                                      const double tmp51 = w2*(A_01_0 - A_10_1);
2554                                      const double tmp32_1 = tmp13_0*w2;                                      const double tmp52 = w5*(-A_01_2 - A_10_1);
2555                                      const double tmp33_1 = tmp14_0*w5;                                      const double tmp53 = w2*(A_01_1 + A_10_2);
2556                                      const double tmp34_1 = tmp15_0*w8;                                      const double tmp54 = w5*(-A_01_2 - A_10_2);
2557                                      const double tmp35_1 = A_01_0*w8;                                      const double tmp55 = w2*(A_01_1 + A_10_1);
2558                                      const double tmp36_1 = tmp16_0*w1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp9;
2559                                      const double tmp37_1 = A_10_1*w4;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;
2560                                      const double tmp38_1 = tmp5_0*w15;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp31 + tmp34 + tmp38 + tmp39 + tmp40 + tmp41;
2561                                      const double tmp39_1 = A_10_2*w7;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp28 + tmp52 + tmp53 + tmp7 + tmp8;
2562                                      const double tmp40_1 = tmp11_0*w17;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp2 + tmp3 + tmp31 + tmp50 + tmp51;
2563                                      const double tmp41_1 = A_01_3*w2;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp16 + tmp17 + tmp21 + tmp25 + tmp28 + tmp54 + tmp55;
2564                                      const double tmp42_1 = tmp12_0*w16;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10 + tmp6 + tmp7 + tmp8 + tmp9;
2565                                      const double tmp43_1 = tmp17_0*w5;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1 + tmp30 + tmp33 + tmp34 + tmp44 + tmp45;
2566                                      const double tmp44_1 = tmp7_0*w7;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1 + tmp34 + tmp38 + tmp40 + tmp42 + tmp43;
2567                                      const double tmp45_1 = tmp6_0*w4;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp36 + tmp37 + tmp7 + tmp8 + tmp9;
2568                                      const double tmp46_1 = A_01_3*w8;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp15 + tmp19 + tmp22 + tmp23 + tmp28 + tmp48 + tmp49;
2569                                      const double tmp47_1 = A_10_2*w4;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0 + tmp11 + tmp12 + tmp31 + tmp46 + tmp47;
2570                                      const double tmp48_1 = A_10_1*w7;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp27 + tmp28 + tmp29 + tmp7 + tmp8;
2571                                      const double tmp49_1 = tmp12_0*w17;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35;
2572                                      const double tmp50_1 = A_01_0*w2;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0 + tmp1 + tmp11 + tmp12 + tmp13 + tmp14;
2573                                      const double tmp51_1 = tmp11_0*w16;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp9;
                                     const double tmp52_1 = tmp18_0*w8;  
                                     const double tmp53_1 = tmp19_0*w2;  
                                     const double tmp54_1 = tmp13_0*w8;  
                                     const double tmp55_1 = tmp15_0*w2;  
                                     const double tmp56_1 = tmp19_0*w8;  
                                     const double tmp57_1 = tmp18_0*w2;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
2574                                  }                                  }
2575                              }                              }
2576                          } else { // constant data                          } else { // constant data
# Line 2593  void Rectangle::assemblePDESystem(Paso_S Line 2580  void Rectangle::assemblePDESystem(Paso_S
2580                                      const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];                                      const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
2581                                      const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];                                      const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
2582                                      const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];                                      const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
2583                                      const double tmp0_0 = A_01 + A_10;                                      const double tmp0 = w9*(-A_01 - A_10);
2584                                      const double tmp0_1 = A_00*w18;                                      const double tmp1 = w9*(-A_01 + A_10);
2585                                      const double tmp1_1 = A_01*w19;                                      const double tmp2 = w9*(A_01 + A_10);
2586                                      const double tmp2_1 = A_10*w20;                                      const double tmp3 = w9*(A_01 - A_10);
2587                                      const double tmp3_1 = A_11*w21;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;
2588                                      const double tmp4_1 = A_00*w22;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*A_00*w6 + 6*A_01*w1 + A_10*w9 + 4*A_11*w3;
2589                                      const double tmp5_1 = tmp0_0*w19;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;
2590                                      const double tmp6_1 = A_11*w23;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp2;
2591                                      const double tmp7_1 = A_11*w25;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp3;
2592                                      const double tmp8_1 = A_00*w24;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;
2593                                      const double tmp9_1 = tmp0_0*w20;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp0;
2594                                      const double tmp10_1 = A_01*w20;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;
2595                                      const double tmp11_1 = A_11*w27;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;
2596                                      const double tmp12_1 = A_00*w26;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp0;
2597                                      const double tmp13_1 = A_10*w19;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;
2598                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp3;
2599                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp2;
2600                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;
2601                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp1;
2602                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
2603                                  }                                  }
2604                              }                              }
2605                          }                          }
# Line 2645  void Rectangle::assemblePDESystem(Paso_S Line 2621  void Rectangle::assemblePDESystem(Paso_S
2621                                      const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];                                      const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
2622                                      const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];                                      const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
2623                                      const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];                                      const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
2624                                      const double tmp0_0 = B_1_0 + B_1_1;                                      const double tmp0 = w15*(B_1_2 + B_1_3);
2625                                      const double tmp1_0 = B_1_2 + B_1_3;                                      const double tmp1 = w12*(B_1_0 + B_1_1);
2626                                      const double tmp2_0 = B_0_1 + B_0_3;                                      const double tmp2 = w15*(B_1_0 + B_1_1);
2627                                      const double tmp3_0 = B_0_0 + B_0_2;                                      const double tmp3 = w16*(-B_0_1 - B_0_3);
2628                                      const double tmp63_1 = B_1_1*w42;                                      const double tmp4 = w11*(-B_0_0 - B_0_2);
2629                                      const double tmp79_1 = B_1_1*w40;                                      const double tmp5 = w12*(B_1_2 + B_1_3);
2630                                      const double tmp37_1 = tmp3_0*w35;                                      const double tmp6 = w15*(-B_1_0 - B_1_1);
2631                                      const double tmp8_1 = tmp0_0*w32;                                      const double tmp7 = w12*(-B_1_2 - B_1_3);
2632                                      const double tmp71_1 = B_0_1*w34;                                      const double tmp8 = w15*(-B_1_2 - B_1_3);
2633                                      const double tmp19_1 = B_0_3*w31;                                      const double tmp9 = w12*(-B_1_0 - B_1_1);
2634                                      const double tmp15_1 = B_0_3*w34;                                      const double tmp10 = w11*(-B_0_1 - B_0_3);
2635                                      const double tmp9_1 = tmp3_0*w34;                                      const double tmp11 = w16*(-B_0_0 - B_0_2);
2636                                      const double tmp35_1 = B_1_0*w36;                                      const double tmp12 = w16*(B_0_0 + B_0_2);
2637                                      const double tmp66_1 = B_0_3*w28;                                      const double tmp13 = w11*(B_0_1 + B_0_3);
2638                                      const double tmp28_1 = B_1_0*w42;                                      const double tmp14 = w11*(B_0_0 + B_0_2);
2639                                      const double tmp22_1 = B_1_0*w40;                                      const double tmp15 = w16*(B_0_1 + B_0_3);
2640                                      const double tmp16_1 = B_1_2*w29;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=B_0_0*w13 + B_0_1*w11 + B_0_2*w16 + B_0_3*w14 + B_1_0*w17 + B_1_1*w15 + B_1_2*w12 + B_1_3*w18;
2641                                      const double tmp6_1 = tmp2_0*w35;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=B_0_0*w11 + B_0_1*w13 + B_0_2*w14 + B_0_3*w16 + tmp0 + tmp1;
2642                                      const double tmp55_1 = B_1_3*w40;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=B_1_0*w12 + B_1_1*w18 + B_1_2*w17 + B_1_3*w15 + tmp14 + tmp15;
2643                                      const double tmp50_1 = B_1_3*w42;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp2 + tmp5;
2644                                      const double tmp7_1 = tmp1_0*w29;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-B_0_0*w13 - B_0_1*w11 - B_0_2*w16 - B_0_3*w14 + tmp0 + tmp1;
2645                                      const double tmp1_1 = tmp1_0*w32;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-B_0_0*w11 - B_0_1*w13 - B_0_2*w14 - B_0_3*w16 + B_1_0*w15 + B_1_1*w17 + B_1_2*w18 + B_1_3*w12;
2646                                      const double tmp57_1 = B_0_3*w30;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2647                                      const double tmp18_1 = B_1_1*w32;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=B_1_0*w18 + B_1_1*w12 + B_1_2*w15 + B_1_3*w17 + tmp10 + tmp11;
2648                                      const double tmp53_1 = B_1_0*w41;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-B_1_0*w17 - B_1_1*w15 - B_1_2*w12 - B_1_3*w18 + tmp14 + tmp15;
2649                                      const double tmp61_1 = B_1_3*w36;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2650                                      const double tmp27_1 = B_0_3*w38;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=B_0_0*w16 + B_0_1*w14 + B_0_2*w13 + B_0_3*w11 - B_1_0*w12 - B_1_1*w18 - B_1_2*w17 - B_1_3*w15;
2651                                      const double tmp64_1 = B_0_2*w30;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=B_0_0*w14 + B_0_1*w16 + B_0_2*w11 + B_0_3*w13 + tmp6 + tmp7;
2652                                      const double tmp76_1 = B_0_1*w38;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp3 + tmp4 + tmp8 + tmp9;
2653                                      const double tmp39_1 = tmp2_0*w34;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-B_1_0*w15 - B_1_1*w17 - B_1_2*w18 - B_1_3*w12 + tmp10 + tmp11;
2654                                      const double tmp62_1 = B_0_1*w31;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-B_0_0*w16 - B_0_1*w14 - B_0_2*w13 - B_0_3*w11 + tmp6 + tmp7;
2655                                      const double tmp56_1 = B_0_0*w31;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-B_0_0*w14 - B_0_1*w16 - B_0_2*w11 - B_0_3*w13 - B_1_0*w18 - B_1_1*w12 - B_1_2*w15 - B_1_3*w17;
                                     const double tmp49_1 = B_1_1*w36;  
                                     const double tmp2_1 = B_0_2*w31;  
                                     const double tmp23_1 = B_0_2*w33;  
                                     const double tmp38_1 = B_1_1*w43;  
                                     const double tmp74_1 = B_1_2*w41;  
                                     const double tmp43_1 = B_1_1*w41;  
                                     const double tmp58_1 = B_0_2*w28;  
                                     const double tmp67_1 = B_0_0*w33;  
                                     const double tmp33_1 = tmp0_0*w39;  
                                     const double tmp4_1 = B_0_0*w28;  
                                     const double tmp20_1 = B_0_0*w30;  
                                     const double tmp13_1 = B_0_2*w38;  
                                     const double tmp65_1 = B_1_2*w43;  
                                     const double tmp0_1 = tmp0_0*w29;  
                                     const double tmp41_1 = tmp3_0*w33;  
                                     const double tmp73_1 = B_0_2*w37;  
                                     const double tmp69_1 = B_0_0*w38;  
                                     const double tmp48_1 = B_1_2*w39;  
                                     const double tmp59_1 = B_0_1*w33;  
                                     const double tmp17_1 = B_1_3*w41;  
                                     const double tmp5_1 = B_0_3*w33;  
                                     const double tmp3_1 = B_0_1*w30;  
                                     const double tmp21_1 = B_0_1*w28;  
                                     const double tmp42_1 = B_1_0*w29;  
                                     const double tmp54_1 = B_1_2*w32;  
                                     const double tmp60_1 = B_1_0*w39;  
                                     const double tmp32_1 = tmp1_0*w36;  
                                     const double tmp10_1 = B_0_1*w37;  
                                     const double tmp14_1 = B_0_0*w35;  
                                     const double tmp29_1 = B_0_1*w35;  
                                     const double tmp26_1 = B_1_2*w36;  
                                     const double tmp30_1 = B_1_3*w43;  
                                     const double tmp70_1 = B_0_2*w35;  
                                     const double tmp34_1 = B_1_3*w39;  
                                     const double tmp51_1 = B_1_0*w43;  
                                     const double tmp31_1 = B_0_2*w34;  
                                     const double tmp45_1 = tmp3_0*w28;  
                                     const double tmp11_1 = tmp1_0*w39;  
                      &n