/[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 3766 by caltinay, Thu Jan 12 08:17:49 2012 UTC trunk/ripley/src/Rectangle.cpp revision 4013 by caltinay, Thu Oct 4 03:13:27 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2012 by University of Queensland  * Copyright (c) 2003-2012 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>
17  extern "C" {  extern "C" {
18  #include "paso/SystemMatrix.h"  #include <paso/SystemMatrix.h>
19  }  }
20    
21    #ifdef USE_NETCDF
22    #include <netcdfcpp.h>
23    #endif
24    
25  #if USE_SILO  #if USE_SILO
26  #include <silo.h>  #include <silo.h>
27  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 29  using namespace std; Line 35  using namespace std;
35    
36  namespace ripley {  namespace ripley {
37    
38  Rectangle::Rectangle(int n0, int n1, double l0, double l1, int d0, int d1) :  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
39                         double y1, int d0, int d1) :
40      RipleyDomain(2),      RipleyDomain(2),
41      m_gNE0(n0),      m_x0(x0),
42      m_gNE1(n1),      m_y0(y0),
43      m_l0(l0),      m_l0(x1-x0),
44      m_l1(l1),      m_l1(y1-y0)
     m_NX(d0),  
     m_NY(d1)  
45  {  {
46        // ignore subdivision parameters for serial run
47        if (m_mpiInfo->size == 1) {
48            d0=1;
49            d1=1;
50        }
51    
52        bool warn=false;
53        // if number of subdivisions is non-positive, try to subdivide by the same
54        // ratio as the number of elements
55        if (d0<=0 && d1<=0) {
56            warn=true;
57            d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
58            d1=m_mpiInfo->size/d0;
59            if (d0*d1 != m_mpiInfo->size) {
60                // ratios not the same so subdivide side with more elements only
61                if (n0>n1) {
62                    d0=0;
63                    d1=1;
64                } else {
65                    d0=1;
66                    d1=0;
67                }
68            }
69        }
70        if (d0<=0) {
71            // d1 is preset, determine d0 - throw further down if result is no good
72            d0=m_mpiInfo->size/d1;
73        } else if (d1<=0) {
74            // d0 is preset, determine d1 - throw further down if result is no good
75            d1=m_mpiInfo->size/d0;
76        }
77    
78        m_NX=d0;
79        m_NY=d1;
80    
81      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
82      // among number of ranks      // among number of ranks
83      if (m_NX*m_NY != m_mpiInfo->size)      if (m_NX*m_NY != m_mpiInfo->size)
84          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
85    
86      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)      if (warn) {
87          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
88                << d1 << "). This may not be optimal!" << endl;
89        }
90    
91        if ((n0+1)%m_NX > 0) {
92            double Dx=m_l0/n0;
93            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
94            m_l0=Dx*n0;
95            cout << "Warning: Adjusted number of elements and length. N0="
96                << n0 << ", l0=" << m_l0 << endl;
97        }
98        if ((n1+1)%m_NY > 0) {
99            double Dy=m_l1/n1;
100            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
101            m_l1=Dy*n1;
102            cout << "Warning: Adjusted number of elements and length. N1="
103                << n1 << ", l1=" << m_l1 << endl;
104        }
105    
106        m_gNE0=n0;
107        m_gNE1=n1;
108    
109      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
110          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
# Line 76  Rectangle::Rectangle(int n0, int n1, dou Line 136  Rectangle::Rectangle(int n0, int n1, dou
136    
137      populateSampleIds();      populateSampleIds();
138      createPattern();      createPattern();
   
     if (m_gNE0>2*m_NX && m_gNE1>2*m_NY && m_gNE0%2==0 && m_gNE1%2==0) {  
         m_coarseMesh.reset(new Rectangle(n0/2, n1/2, l0, l1, d0, d1));  
     }  
139  }  }
140    
141  Rectangle::~Rectangle()  Rectangle::~Rectangle()
142  {  {
143        Paso_SystemMatrixPattern_free(m_pattern);
144        Paso_Connector_free(m_connector);
145  }  }
146    
147  string Rectangle::getDescription() const  string Rectangle::getDescription() const
# Line 97  bool Rectangle::operator==(const Abstrac Line 155  bool Rectangle::operator==(const Abstrac
155      if (o) {      if (o) {
156          return (RipleyDomain::operator==(other) &&          return (RipleyDomain::operator==(other) &&
157                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1
158                    && m_x0==o->m_x0 && m_y0==o->m_y0
159                  && m_l0==o->m_l0 && m_l1==o->m_l1                  && m_l0==o->m_l0 && m_l1==o->m_l1
160                  && m_NX==o->m_NX && m_NY==o->m_NY);                  && m_NX==o->m_NX && m_NY==o->m_NY);
161      }      }
# Line 104  bool Rectangle::operator==(const Abstrac Line 163  bool Rectangle::operator==(const Abstrac
163      return false;      return false;
164  }  }
165    
166    void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
167                const vector<int>& first, const vector<int>& numValues) const
168    {
169    #ifdef USE_NETCDF
170        // check destination function space
171        int myN0, myN1;
172        if (out.getFunctionSpace().getTypeCode() == Nodes) {
173            myN0 = m_N0;
174            myN1 = m_N1;
175        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
176                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
177            myN0 = m_NE0;
178            myN1 = m_NE1;
179        } else
180            throw RipleyException("readNcGrid(): invalid function space for output data object");
181    
182        if (first.size() != 2)
183            throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
184    
185        if (numValues.size() != 2)
186            throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
187    
188        // check file existence and size
189        NcFile f(filename.c_str(), NcFile::ReadOnly);
190        if (!f.is_valid())
191            throw RipleyException("readNcGrid(): cannot open file");
192    
193        NcVar* var = f.get_var(varname.c_str());
194        if (!var)
195            throw RipleyException("readNcGrid(): invalid variable");
196    
197        // TODO: rank>0 data support
198        const int numComp = out.getDataPointSize();
199        if (numComp > 1)
200            throw RipleyException("readNcGrid(): only scalar data supported");
201    
202        const int dims = var->num_dims();
203        const long *edges = var->edges();
204    
205        // is this a slice of the data object (dims!=2)?
206        // note the expected ordering of edges (as in numpy: y,x)
207        if ( (dims==2 && (numValues[1] > edges[0] || numValues[0] > edges[1]))
208                || (dims==1 && numValues[1]>1) ) {
209            throw RipleyException("readNcGrid(): not enough data in file");
210        }
211    
212        // check if this rank contributes anything
213        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
214                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1)
215            return;
216    
217        // now determine how much this rank has to write
218    
219        // first coordinates in data object to write to
220        const int first0 = max(0, first[0]-m_offset0);
221        const int first1 = max(0, first[1]-m_offset1);
222        // indices to first value in file
223        const int idx0 = max(0, m_offset0-first[0]);
224        const int idx1 = max(0, m_offset1-first[1]);
225        // number of values to write
226        const int num0 = min(numValues[0]-idx0, myN0-first0);
227        const int num1 = min(numValues[1]-idx1, myN1-first1);
228    
229        vector<double> values(num0*num1);
230        if (dims==2) {
231            var->set_cur(idx1, idx0);
232            var->get(&values[0], num1, num0);
233        } else {
234            var->set_cur(idx0);
235            var->get(&values[0], num0);
236        }
237    
238        const int dpp = out.getNumDataPointsPerSample();
239        out.requireWrite();
240    
241        for (index_t y=0; y<num1; y++) {
242    #pragma omp parallel for
243            for (index_t x=0; x<num0; x++) {
244                const int dataIndex = (first1+y)*myN0+first0+x;
245                const int srcIndex=y*num0+x;
246                double* dest = out.getSampleDataRW(dataIndex);
247                for (index_t q=0; q<dpp; q++) {
248                    *dest++ = values[srcIndex];
249                }
250            }
251        }
252    #else
253        throw RipleyException("readNcGrid(): not compiled with netCDF support");
254    #endif
255    }
256    
257    void Rectangle::readBinaryGrid(escript::Data& out, string filename,
258                const vector<int>& first, const vector<int>& numValues) const
259    {
260        // check destination function space
261        int myN0, myN1;
262        if (out.getFunctionSpace().getTypeCode() == Nodes) {
263            myN0 = m_N0;
264            myN1 = m_N1;
265        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
266                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
267            myN0 = m_NE0;
268            myN1 = m_NE1;
269        } else
270            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
271    
272        // check file existence and size
273        ifstream f(filename.c_str(), ifstream::binary);
274        if (f.fail()) {
275            throw RipleyException("readBinaryGrid(): cannot open file");
276        }
277        f.seekg(0, ios::end);
278        const int numComp = out.getDataPointSize();
279        const int filesize = f.tellg();
280        const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(float);
281        if (filesize < reqsize) {
282            f.close();
283            throw RipleyException("readBinaryGrid(): not enough data in file");
284        }
285    
286        // check if this rank contributes anything
287        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
288                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1) {
289            f.close();
290            return;
291        }
292    
293        // now determine how much this rank has to write
294    
295        // first coordinates in data object to write to
296        const int first0 = max(0, first[0]-m_offset0);
297        const int first1 = max(0, first[1]-m_offset1);
298        // indices to first value in file
299        const int idx0 = max(0, m_offset0-first[0]);
300        const int idx1 = max(0, m_offset1-first[1]);
301        // number of values to write
302        const int num0 = min(numValues[0]-idx0, myN0-first0);
303        const int num1 = min(numValues[1]-idx1, myN1-first1);
304    
305        out.requireWrite();
306        vector<float> values(num0*numComp);
307        const int dpp = out.getNumDataPointsPerSample();
308    
309        for (index_t y=0; y<num1; y++) {
310            const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);
311            f.seekg(fileofs*sizeof(float));
312            f.read((char*)&values[0], num0*numComp*sizeof(float));
313            for (index_t x=0; x<num0; x++) {
314                double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0);
315                for (index_t c=0; c<numComp; c++) {
316                    for (index_t q=0; q<dpp; q++) {
317                        *dest++ = static_cast<double>(values[x*numComp+c]);
318                    }
319                }
320            }
321        }
322    
323        f.close();
324    }
325    
326  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
327  {  {
328  #if USE_SILO  #if USE_SILO
# Line 112  void Rectangle::dump(const string& fileN Line 331  void Rectangle::dump(const string& fileN
331          fn+=".silo";          fn+=".silo";
332      }      }
333    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
334      int driver=DB_HDF5;          int driver=DB_HDF5;    
335      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
336        const char* blockDirFmt = "/block%04d";
337    
338  #ifdef ESYS_MPI  #ifdef ESYS_MPI
339      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
340        const int NUM_SILO_FILES = 1;
341  #endif  #endif
342    
343      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 244  void Rectangle::dump(const string& fileN Line 463  void Rectangle::dump(const string& fileN
463      }      }
464    
465  #else // USE_SILO  #else // USE_SILO
466      throw RipleyException("dump(): no Silo support");      throw RipleyException("dump: no Silo support");
467  #endif  #endif
468  }  }
469    
# Line 252  const int* Rectangle::borrowSampleRefere Line 471  const int* Rectangle::borrowSampleRefere
471  {  {
472      switch (fsType) {      switch (fsType) {
473          case Nodes:          case Nodes:
474            case ReducedNodes: // FIXME: reduced
475              return &m_nodeId[0];              return &m_nodeId[0];
         case ReducedNodes:  
             if (m_coarseMesh)  
                 return m_coarseMesh->borrowSampleReferenceIDs(Nodes);  
             break;  
476          case DegreesOfFreedom:          case DegreesOfFreedom:
477            case ReducedDegreesOfFreedom: // FIXME: reduced
478              return &m_dofId[0];              return &m_dofId[0];
         case ReducedDegreesOfFreedom:  
             if (m_coarseMesh)  
                 return m_coarseMesh->borrowSampleReferenceIDs(DegreesOfFreedom);  
             break;  
479          case Elements:          case Elements:
480          case ReducedElements:          case ReducedElements:
481              return &m_elementId[0];              return &m_elementId[0];
# Line 274  const int* Rectangle::borrowSampleRefere Line 487  const int* Rectangle::borrowSampleRefere
487      }      }
488    
489      stringstream msg;      stringstream msg;
490      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
491      throw RipleyException(msg.str());      throw RipleyException(msg.str());
492  }  }
493    
# Line 286  bool Rectangle::ownSample(int fsType, in Line 498  bool Rectangle::ownSample(int fsType, in
498    
499      switch (fsType) {      switch (fsType) {
500          case Nodes:          case Nodes:
501            case ReducedNodes: // FIXME: reduced
502              return (m_dofMap[id] < getNumDOF());              return (m_dofMap[id] < getNumDOF());
         case ReducedNodes:  
             if (m_coarseMesh)  
                 return m_coarseMesh->ownSample(Nodes, id);  
             break;  
503          case DegreesOfFreedom:          case DegreesOfFreedom:
504          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
505              return true;              return true;
# Line 302  bool Rectangle::ownSample(int fsType, in Line 511  bool Rectangle::ownSample(int fsType, in
511          case ReducedFaceElements:          case ReducedFaceElements:
512              {              {
513                  // determine which face the sample belongs to before                  // determine which face the sample belongs to before
514                  // checking ownership of face element's first node                  // checking ownership of corresponding element's first node
515                  const IndexVector faces = getNumFacesPerBoundary();                  const IndexVector faces = getNumFacesPerBoundary();
516                  dim_t n=0;                  dim_t n=0;
517                  for (size_t i=0; i<faces.size(); i++) {                  for (size_t i=0; i<faces.size(); i++) {
# Line 310  bool Rectangle::ownSample(int fsType, in Line 519  bool Rectangle::ownSample(int fsType, in
519                      if (id<n) {                      if (id<n) {
520                          index_t k;                          index_t k;
521                          if (i==1)                          if (i==1)
522                              k=m_N0-1;                              k=m_N0-2;
523                          else if (i==3)                          else if (i==3)
524                              k=m_N0*(m_N1-1);                              k=m_N0*(m_N1-2);
525                          else                          else
526                              k=0;                              k=0;
527                          // determine whether to move right or up                          // determine whether to move right or up
# Line 327  bool Rectangle::ownSample(int fsType, in Line 536  bool Rectangle::ownSample(int fsType, in
536      }      }
537    
538      stringstream msg;      stringstream msg;
539      msg << "ownSample() not implemented for "      msg << "ownSample: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
540      throw RipleyException(msg.str());      throw RipleyException(msg.str());
541  }  }
542    
# Line 429  void Rectangle::setToNormal(escript::Dat Line 637  void Rectangle::setToNormal(escript::Dat
637    
638      } else {      } else {
639          stringstream msg;          stringstream msg;
640          msg << "setToNormal() not implemented for "          msg << "setToNormal: invalid function space type "
641              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());              << out.getFunctionSpace().getTypeCode();
642          throw RipleyException(msg.str());          throw RipleyException(msg.str());
643      }      }
644  }  }
# Line 443  void Rectangle::setToSize(escript::Data& Line 651  void Rectangle::setToSize(escript::Data&
651          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
652          const double hSize=getFirstCoordAndSpacing(0).second;          const double hSize=getFirstCoordAndSpacing(0).second;
653          const double vSize=getFirstCoordAndSpacing(1).second;          const double vSize=getFirstCoordAndSpacing(1).second;
654          const double size=min(hSize,vSize);          const double size=sqrt(hSize*hSize+vSize*vSize);
655  #pragma omp parallel for  #pragma omp parallel for
656          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
657              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 492  void Rectangle::setToSize(escript::Data& Line 700  void Rectangle::setToSize(escript::Data&
700    
701      } else {      } else {
702          stringstream msg;          stringstream msg;
703          msg << "setToSize() not implemented for "          msg << "setToSize: invalid function space type "
704              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());              << out.getFunctionSpace().getTypeCode();
705          throw RipleyException(msg.str());          throw RipleyException(msg.str());
706      }      }
707  }  }
# Line 501  void Rectangle::setToSize(escript::Data& Line 709  void Rectangle::setToSize(escript::Data&
709  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
710                                                  bool reducedColOrder) const                                                  bool reducedColOrder) const
711  {  {
712        /* FIXME: reduced
713      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
714          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
715        */
716      return m_pattern;      return m_pattern;
717  }  }
718    
# Line 569  IndexVector Rectangle::getNumSubdivision Line 778  IndexVector Rectangle::getNumSubdivision
778  pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
779  {  {
780      if (dim==0) {      if (dim==0) {
781          return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);          return pair<double,double>(m_x0+(m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
782      } else if (dim==1) {      } else if (dim==1) {
783          return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);          return pair<double,double>(m_y0+(m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
784      }      }
785      throw RipleyException("getFirstCoordAndSpacing(): invalid argument");      throw RipleyException("getFirstCoordAndSpacing: invalid argument");
786  }  }
787    
788  //protected  //protected
# Line 640  void Rectangle::assembleGradient(escript Line 849  void Rectangle::assembleGradient(escript
849    
850      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
851          out.requireWrite();          out.requireWrite();
852  #pragma omp parallel for  #pragma omp parallel
853          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
854              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
855                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
856                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
857                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
858                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
859                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
860                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
861                      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_N0)), numComp*sizeof(double));
862                      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_N0)), numComp*sizeof(double));
863                      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_N0)), numComp*sizeof(double));
864                      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_N0)), numComp*sizeof(double));
865                      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_NE0));
866                      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) {
867                      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_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
868                      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_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
869                  } /* end of component loop i */                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
870              } /* end of k0 loop */                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
871          } /* end of k1 loop */                          o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
872                            o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
873                            o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
874                            o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
875                        } // end of component loop i
876                    } // end of k0 loop
877                } // end of k1 loop
878            } // end of parallel section
879      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
880          out.requireWrite();          out.requireWrite();
881  #pragma omp parallel for  #pragma omp parallel
882          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
883              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
884                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
885                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
886                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
887                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
888                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
889                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
890                      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_N0)), numComp*sizeof(double));
891                      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_N0)), numComp*sizeof(double));
892                  } /* end of component loop i */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
893              } /* end of k0 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
894          } /* end of k1 loop */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
895                        for (index_t i=0; i < numComp; ++i) {
896                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
897                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
898                        } // end of component loop i
899                    } // end of k0 loop
900                } // end of k1 loop
901            } // end of parallel section
902      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
903          out.requireWrite();          out.requireWrite();
904  #pragma omp parallel  #pragma omp parallel
905          {          {
906                vector<double> f_00(numComp);
907                vector<double> f_01(numComp);
908                vector<double> f_10(numComp);
909                vector<double> f_11(numComp);
910              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
911  #pragma omp for nowait  #pragma omp for nowait
912                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
913                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
914                      const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
915                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
916                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
917                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
918                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
919                          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_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
920                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
921                          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_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
922                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
923                      } /* end of component loop i */                      } // end of component loop i
924                  } /* end of k1 loop */                  } // end of k1 loop
925              } /* end of face 0 */              } // end of face 0
926              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
927  #pragma omp for nowait  #pragma omp for nowait
928                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
929                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
930                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
931                      const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
932                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
933                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
934                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
935                          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_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
936                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
937                          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_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
938                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
939                      } /* end of component loop i */                      } // end of component loop i
940                  } /* end of k1 loop */                  } // end of k1 loop
941              } /* end of face 1 */              } // end of face 1
942              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
943  #pragma omp for nowait  #pragma omp for nowait
944                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
945                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
946                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
947                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
948                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
949                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
950                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
951                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
952                          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_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
953                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
954                          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_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
955                      } /* end of component loop i */                      } // end of component loop i
956                  } /* end of k0 loop */                  } // end of k0 loop
957              } /* end of face 2 */              } // end of face 2
958              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
959  #pragma omp for nowait  #pragma omp for nowait
960                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
961                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
962                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
963                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
964                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
965                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
966                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
967                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
968                          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_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
969                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
970                          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_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
971                      } /* end of component loop i */                      } // end of component loop i
972                  } /* end of k0 loop */                  } // end of k0 loop
973              } /* end of face 3 */              } // end of face 3
974          } // end of parallel section          } // end of parallel section
975    
976      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
977          out.requireWrite();          out.requireWrite();
978  #pragma omp parallel  #pragma omp parallel
979          {          {
980                vector<double> f_00(numComp);
981                vector<double> f_01(numComp);
982                vector<double> f_10(numComp);
983                vector<double> f_11(numComp);
984              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
985  #pragma omp for nowait  #pragma omp for nowait
986                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
987                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
988                      const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
989                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
990                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
991                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
992                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
993                          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)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
994                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
995                      } /* end of component loop i */                      } // end of component loop i
996                  } /* end of k1 loop */                  } // end of k1 loop
997              } /* end of face 0 */              } // end of face 0
998              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
999  #pragma omp for nowait  #pragma omp for nowait
1000                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1001                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
1002                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
1003                      const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
1004                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
1005                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1006                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1007                          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)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
1008                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
1009                      } /* end of component loop i */                      } // end of component loop i
1010                  } /* end of k1 loop */                  } // end of k1 loop
1011              } /* end of face 1 */              } // end of face 1
1012              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1013  #pragma omp for nowait  #pragma omp for nowait
1014                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1015                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
1016                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
1017                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
1018                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
1019                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1020                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1021                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
1022                          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)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
1023                      } /* end of component loop i */                      } // end of component loop i
1024                  } /* end of k0 loop */                  } // end of k0 loop
1025              } /* end of face 2 */              } // end of face 2
1026              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1027  #pragma omp for nowait  #pragma omp for nowait
1028                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1029                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
1030                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
1031                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
1032                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
1033                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1034                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1035                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
1036                          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)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
1037                      } /* end of component loop i */                      } // end of component loop i
1038                  } /* end of k0 loop */                  } // end of k0 loop
1039              } /* end of face 3 */              } // end of face 3
1040          } // end of parallel section          } // end of parallel section
1041      }      }
1042  }  }
# Line 817  void Rectangle::assembleIntegrate(vector Line 1049  void Rectangle::assembleIntegrate(vector
1049      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
1050      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1051      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1052      if (arg.getFunctionSpace().getTypeCode() == Elements) {      const int fs=arg.getFunctionSpace().getTypeCode();
1053          const double w = h0*h1/4.;      if (fs == Elements && arg.actsExpanded()) {
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 = h0*h1/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_ownNE1; ++k1) {
1060                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1061                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
1062                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1063                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1064                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1065                          const register double f2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
1066                          const register 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        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1077          const double w = h0*h1;          const double w = h0*h1;
1078  #pragma omp parallel  #pragma omp parallel
1079          {          {
# Line 851  void Rectangle::assembleIntegrate(vector Line 1084  void Rectangle::assembleIntegrate(vector
1084                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
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 = h0/2.;
1100                const double w1 = h1/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_ownNE1; ++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 register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1107                          const register 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) {
# Line 882  void Rectangle::assembleIntegrate(vector Line 1115  void Rectangle::assembleIntegrate(vector
1115                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++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 register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1119                          const register 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) {
# Line 894  void Rectangle::assembleIntegrate(vector Line 1127  void Rectangle::assembleIntegrate(vector
1127                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++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 register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1131                          const register 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) {
# Line 906  void Rectangle::assembleIntegrate(vector Line 1139  void Rectangle::assembleIntegrate(vector
1139                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++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 register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1143                          const register 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);
# Line 927  void Rectangle::assembleIntegrate(vector Line 1160  void Rectangle::assembleIntegrate(vector
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]*h1;
1163                      }  /* end of component loop i */                      }
1164                  } /* end of k1 loop */                  }
1165              }              }
1166    
1167              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
# Line 937  void Rectangle::assembleIntegrate(vector Line 1170  void Rectangle::assembleIntegrate(vector
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]*h1;
1173                      }  /* end of component loop i */                      }
1174                  } /* end of k1 loop */                  }
1175              }              }
1176    
1177              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
# Line 947  void Rectangle::assembleIntegrate(vector Line 1180  void Rectangle::assembleIntegrate(vector
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]*h0;
1183                      }  /* end of component loop i */                      }
1184                  } /* end of k0 loop */                  }
1185              }              }
1186    
1187              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
# Line 957  void Rectangle::assembleIntegrate(vector Line 1190  void Rectangle::assembleIntegrate(vector
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]*h0;
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
# Line 1035  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
# Line 1285  void Rectangle::createPattern() Line 1519  void Rectangle::createPattern()
1519      Paso_Pattern_free(rowPattern);      Paso_Pattern_free(rowPattern);
1520  }  }
1521    
1522    //private
1523    void Rectangle::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
1524             const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
1525             bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
1526    {
1527        IndexVector rowIndex;
1528        rowIndex.push_back(m_dofMap[firstNode]);
1529        rowIndex.push_back(m_dofMap[firstNode+1]);
1530        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
1531        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
1532        if (addF) {
1533            double *F_p=F.getSampleDataRW(0);
1534            for (index_t i=0; i<rowIndex.size(); i++) {
1535                if (rowIndex[i]<getNumDOF()) {
1536                    for (index_t eq=0; eq<nEq; eq++) {
1537                        F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
1538                    }
1539                }
1540            }
1541        }
1542        if (addS) {
1543            addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
1544        }
1545    }
1546    
1547  //protected  //protected
1548  void Rectangle::interpolateNodesOnElements(escript::Data& out,  void Rectangle::interpolateNodesOnElements(escript::Data& out,
1549                                          escript::Data& in, bool reduced) const                                          escript::Data& in, bool reduced) const
# Line 1292  void Rectangle::interpolateNodesOnElemen Line 1551  void Rectangle::interpolateNodesOnElemen
1551      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1552      if (reduced) {      if (reduced) {
1553          out.requireWrite();          out.requireWrite();
1554          const double c0 = .25;          const double c0 = 0.25;
1555  #pragma omp parallel for  #pragma omp parallel
1556          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1557              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1558                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_01(numComp);
1559                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1560                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              vector<double> f_11(numComp);
1561                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1562                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1563                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1564                      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_N0)), numComp*sizeof(double));
1565                  } /* end of component loop i */                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
1566              } /* end of k0 loop */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
1567          } /* end of k1 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
1568                        double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1569                        for (index_t i=0; i < numComp; ++i) {
1570                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1571                        } /* end of component loop i */
1572                    } /* end of k0 loop */
1573                } /* end of k1 loop */
1574            } /* end of parallel section */
1575      } else {      } else {
1576          out.requireWrite();          out.requireWrite();
1577          const double c0 = .16666666666666666667;          const double c0 = 0.16666666666666666667;
1578          const double c1 = .044658198738520451079;          const double c1 = 0.044658198738520451079;
1579          const double c2 = .62200846792814621559;          const double c2 = 0.62200846792814621559;
1580  #pragma omp parallel for  #pragma omp parallel
1581          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1582              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1583                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
1584                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1585                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
1586                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1587                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1588                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1589                      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_N0)), numComp*sizeof(double));
1590                      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_N0)), numComp*sizeof(double));
1591                      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_N0)), numComp*sizeof(double));
1592                      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_N0)), numComp*sizeof(double));
1593                  } /* end of component loop i */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1594              } /* end of k0 loop */                      for (index_t i=0; i < numComp; ++i) {
1595          } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1596                            o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1597                            o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1598                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1599                        } /* end of component loop i */
1600                    } /* end of k0 loop */
1601                } /* end of k1 loop */
1602            } /* end of parallel section */
1603      }      }
1604  }  }
1605    
# Line 1337  void Rectangle::interpolateNodesOnFaces( Line 1610  void Rectangle::interpolateNodesOnFaces(
1610      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1611      if (reduced) {      if (reduced) {
1612          out.requireWrite();          out.requireWrite();
1613          const double c0 = .5;          const double c0 = 0.5;
1614  #pragma omp parallel  #pragma omp parallel
1615          {          {
1616                vector<double> f_00(numComp);
1617                vector<double> f_01(numComp);
1618                vector<double> f_10(numComp);
1619                vector<double> f_11(numComp);
1620              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1621  #pragma omp for nowait  #pragma omp for nowait
1622                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1623                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
1624                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
1625                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1626                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1627                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
# Line 1354  void Rectangle::interpolateNodesOnFaces( Line 1631  void Rectangle::interpolateNodesOnFaces(
1631              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1632  #pragma omp for nowait  #pragma omp for nowait
1633                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1634                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
1635                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
1636                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1637                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1638                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
# Line 1365  void Rectangle::interpolateNodesOnFaces( Line 1642  void Rectangle::interpolateNodesOnFaces(
1642              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1643  #pragma omp for nowait  #pragma omp for nowait
1644                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1645                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
1646                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
1647                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1648                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1649                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
# Line 1376  void Rectangle::interpolateNodesOnFaces( Line 1653  void Rectangle::interpolateNodesOnFaces(
1653              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1654  #pragma omp for nowait  #pragma omp for nowait
1655                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1656                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
1657                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
1658                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1659                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1660                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1661                      } /* end of component loop i */                      } /* end of component loop i */
1662                  } /* end of k0 loop */                  } /* end of k0 loop */
1663              } /* end of face 3 */              } /* end of face 3 */
1664          } // end of parallel section          } /* end of parallel section */
1665      } else {      } else {
1666          out.requireWrite();          out.requireWrite();
1667          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1668          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1669  #pragma omp parallel  #pragma omp parallel
1670          {          {
1671                vector<double> f_00(numComp);
1672                vector<double> f_01(numComp);
1673                vector<double> f_10(numComp);
1674                vector<double> f_11(numComp);
1675              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1676  #pragma omp for nowait  #pragma omp for nowait
1677                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1678                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
1679                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
1680                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1681                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1682                          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];
1683                          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];
1684                      } /* end of component loop i */                      } /* end of component loop i */
1685                  } /* end of k1 loop */                  } /* end of k1 loop */
1686              } /* end of face 0 */              } /* end of face 0 */
1687              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1688  #pragma omp for nowait  #pragma omp for nowait
1689                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1690                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
1691                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
1692                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1693                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1694                          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];
1695                          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];
1696                      } /* end of component loop i */                      } /* end of component loop i */
1697                  } /* end of k1 loop */                  } /* end of k1 loop */
1698              } /* end of face 1 */              } /* end of face 1 */
1699              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1700  #pragma omp for nowait  #pragma omp for nowait
1701                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1702                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
1703                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
1704                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1705                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1706                          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];
1707                          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];
1708                      } /* end of component loop i */                      } /* end of component loop i */
1709                  } /* end of k0 loop */                  } /* end of k0 loop */
1710              } /* end of face 2 */              } /* end of face 2 */
1711              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1712  #pragma omp for nowait  #pragma omp for nowait
1713                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1714                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
1715                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
1716                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1717                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1718                          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];
1719                          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];
1720                      } /* end of component loop i */                      } /* end of component loop i */
1721                  } /* end of k0 loop */                  } /* end of k0 loop */
1722              } /* end of face 3 */              } /* end of face 3 */
1723          } // end of parallel section          } /* end of parallel section */
1724      }      }
1725  }  }
1726    
# Line 1447  void Rectangle::interpolateNodesOnFaces( Line 1728  void Rectangle::interpolateNodesOnFaces(
1728  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1729          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1730          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
1731          const escript::Data& X, const escript::Data& Y,          const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
1732  {  {
1733      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1734      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 1471  void Rectangle::assemblePDESingle(Paso_S Line 1751  void Rectangle::assemblePDESingle(Paso_S
1751      const double w16 = -0.01116454968463011277*h0/h1;      const double w16 = -0.01116454968463011277*h0/h1;
1752      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*h0/h1;
1753      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*h1/h0;
1754      const double w19 = 0.25000000000000000000;      const double w19 = 0.25;
1755      const double w20 = -0.25000000000000000000;      const double w20 = -0.25;
1756      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*h0/h1;
1757      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*h1/h0;
1758      const double w23 = -0.16666666666666666667*h0/h1;      const double w23 = -0.16666666666666666667*h0/h1;
# Line 1548  void Rectangle::assemblePDESingle(Paso_S Line 1828  void Rectangle::assemblePDESingle(Paso_S
1828                          add_EM_S=true;                          add_EM_S=true;
1829                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1830                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
1831                              const register 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)];
1832                              const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];                              const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1833                              const register 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)];
1834                              const register 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)];
1835                              const register 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)];
1836                              const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];                              const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1837                              const register 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)];
1838                              const register 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)];
1839                              const register 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)];
1840                              const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];                              const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1841                              const register 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)];
1842                              const register 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)];
1843                              const register 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)];
1844                              const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];                              const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1845                              const register 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)];
1846                              const register 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)];
1847                              const register double tmp0_0 = A_01_0 + A_01_3;                              const double tmp0_0 = A_01_0 + A_01_3;
1848                              const register double tmp1_0 = A_00_0 + A_00_1;                              const double tmp1_0 = A_00_0 + A_00_1;
1849                              const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                              const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1850                              const register double tmp3_0 = A_00_2 + A_00_3;                              const double tmp3_0 = A_00_2 + A_00_3;
1851                              const register double tmp4_0 = A_10_1 + A_10_2;                              const double tmp4_0 = A_10_1 + A_10_2;
1852                              const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                              const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1853                              const register double tmp6_0 = A_01_3 + A_10_0;                              const double tmp6_0 = A_01_3 + A_10_0;
1854                              const register double tmp7_0 = A_01_0 + A_10_3;                              const double tmp7_0 = A_01_0 + A_10_3;
1855                              const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                              const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1856                              const register double tmp9_0 = A_01_0 + A_10_0;                              const double tmp9_0 = A_01_0 + A_10_0;
1857                              const register double tmp12_0 = A_11_0 + A_11_2;                              const double tmp12_0 = A_11_0 + A_11_2;
1858                              const register double tmp10_0 = A_01_3 + A_10_3;                              const double tmp10_0 = A_01_3 + A_10_3;
1859                              const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                              const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
1860                              const register double tmp13_0 = A_01_2 + A_10_1;                              const double tmp13_0 = A_01_2 + A_10_1;
1861                              const register double tmp11_0 = A_11_1 + A_11_3;                              const double tmp11_0 = A_11_1 + A_11_3;
1862                              const register double tmp18_0 = A_01_1 + A_10_1;                              const double tmp18_0 = A_01_1 + A_10_1;
1863                              const register double tmp15_0 = A_01_1 + A_10_2;                              const double tmp15_0 = A_01_1 + A_10_2;
1864                              const register double tmp16_0 = A_10_0 + A_10_3;                              const double tmp16_0 = A_10_0 + A_10_3;
1865                              const register double tmp17_0 = A_01_1 + A_01_2;                              const double tmp17_0 = A_01_1 + A_01_2;
1866                              const register double tmp19_0 = A_01_2 + A_10_2;                              const double tmp19_0 = A_01_2 + A_10_2;
1867                              const register double tmp0_1 = A_10_3*w8;                              const double tmp0_1 = A_10_3*w8;
1868                              const register double tmp1_1 = tmp0_0*w1;                              const double tmp1_1 = tmp0_0*w1;
1869                              const register double tmp2_1 = A_01_1*w4;                              const double tmp2_1 = A_01_1*w4;
1870                              const register double tmp3_1 = tmp1_0*w0;                              const double tmp3_1 = tmp1_0*w0;
1871                              const register double tmp4_1 = A_01_2*w7;                              const double tmp4_1 = A_01_2*w7;
1872                              const register double tmp5_1 = tmp2_0*w3;                              const double tmp5_1 = tmp2_0*w3;
1873                              const register double tmp6_1 = tmp3_0*w6;                              const double tmp6_1 = tmp3_0*w6;
1874                              const register double tmp7_1 = A_10_0*w2;                              const double tmp7_1 = A_10_0*w2;
1875                              const register double tmp8_1 = tmp4_0*w5;                              const double tmp8_1 = tmp4_0*w5;
1876                              const register double tmp9_1 = tmp2_0*w10;                              const double tmp9_1 = tmp2_0*w10;
1877                              const register double tmp14_1 = A_10_0*w8;                              const double tmp14_1 = A_10_0*w8;
1878                              const register double tmp23_1 = tmp3_0*w14;                              const double tmp23_1 = tmp3_0*w14;
1879                              const register double tmp35_1 = A_01_0*w8;                              const double tmp35_1 = A_01_0*w8;
1880                              const register double tmp54_1 = tmp13_0*w8;                              const double tmp54_1 = tmp13_0*w8;
1881                              const register double tmp20_1 = tmp9_0*w4;                              const double tmp20_1 = tmp9_0*w4;
1882                              const register double tmp25_1 = tmp12_0*w12;                              const double tmp25_1 = tmp12_0*w12;
1883                              const register double tmp44_1 = tmp7_0*w7;                              const double tmp44_1 = tmp7_0*w7;
1884                              const register double tmp26_1 = tmp10_0*w4;                              const double tmp26_1 = tmp10_0*w4;
1885                              const register double tmp52_1 = tmp18_0*w8;                              const double tmp52_1 = tmp18_0*w8;
1886                              const register double tmp48_1 = A_10_1*w7;                              const double tmp48_1 = A_10_1*w7;
1887                              const register double tmp46_1 = A_01_3*w8;                              const double tmp46_1 = A_01_3*w8;
1888                              const register double tmp50_1 = A_01_0*w2;                              const double tmp50_1 = A_01_0*w2;
1889                              const register double tmp56_1 = tmp19_0*w8;                              const double tmp56_1 = tmp19_0*w8;
1890                              const register double tmp19_1 = A_10_3*w2;                              const double tmp19_1 = A_10_3*w2;
1891                              const register double tmp47_1 = A_10_2*w4;                              const double tmp47_1 = A_10_2*w4;
1892                              const register double tmp16_1 = tmp3_0*w0;                              const double tmp16_1 = tmp3_0*w0;
1893                              const register double tmp18_1 = tmp1_0*w6;                              const double tmp18_1 = tmp1_0*w6;
1894                              const register double tmp31_1 = tmp11_0*w12;                              const double tmp31_1 = tmp11_0*w12;
1895                              const register double tmp55_1 = tmp15_0*w2;                              const double tmp55_1 = tmp15_0*w2;
1896                              const register double tmp39_1 = A_10_2*w7;                              const double tmp39_1 = A_10_2*w7;
1897                              const register double tmp11_1 = tmp6_0*w7;                              const double tmp11_1 = tmp6_0*w7;
1898                              const register double tmp40_1 = tmp11_0*w17;                              const double tmp40_1 = tmp11_0*w17;
1899                              const register double tmp34_1 = tmp15_0*w8;                              const double tmp34_1 = tmp15_0*w8;
1900                              const register double tmp33_1 = tmp14_0*w5;                              const double tmp33_1 = tmp14_0*w5;
1901                              const register double tmp24_1 = tmp11_0*w13;                              const double tmp24_1 = tmp11_0*w13;
1902                              const register double tmp43_1 = tmp17_0*w5;                              const double tmp43_1 = tmp17_0*w5;
1903                              const register double tmp15_1 = A_01_2*w4;                              const double tmp15_1 = A_01_2*w4;
1904                              const register double tmp53_1 = tmp19_0*w2;                              const double tmp53_1 = tmp19_0*w2;
1905                              const register double tmp27_1 = tmp3_0*w11;                              const double tmp27_1 = tmp3_0*w11;
1906                              const register double tmp32_1 = tmp13_0*w2;                              const double tmp32_1 = tmp13_0*w2;
1907                              const register double tmp10_1 = tmp5_0*w9;                              const double tmp10_1 = tmp5_0*w9;
1908                              const register double tmp37_1 = A_10_1*w4;                              const double tmp37_1 = A_10_1*w4;
1909                              const register double tmp38_1 = tmp5_0*w15;                              const double tmp38_1 = tmp5_0*w15;
1910                              const register double tmp17_1 = A_01_1*w7;                              const double tmp17_1 = A_01_1*w7;
1911                              const register double tmp12_1 = tmp7_0*w4;                              const double tmp12_1 = tmp7_0*w4;
1912                              const register double tmp22_1 = tmp10_0*w7;                              const double tmp22_1 = tmp10_0*w7;
1913                              const register double tmp57_1 = tmp18_0*w2;                              const double tmp57_1 = tmp18_0*w2;
1914                              const register double tmp28_1 = tmp9_0*w7;                              const double tmp28_1 = tmp9_0*w7;
1915                              const register double tmp29_1 = tmp1_0*w14;                              const double tmp29_1 = tmp1_0*w14;
1916                              const register double tmp51_1 = tmp11_0*w16;                              const double tmp51_1 = tmp11_0*w16;
1917                              const register double tmp42_1 = tmp12_0*w16;                              const double tmp42_1 = tmp12_0*w16;
1918                              const register double tmp49_1 = tmp12_0*w17;                              const double tmp49_1 = tmp12_0*w17;
1919                              const register double tmp21_1 = tmp1_0*w11;                              const double tmp21_1 = tmp1_0*w11;
1920                              const register double tmp45_1 = tmp6_0*w4;                              const double tmp45_1 = tmp6_0*w4;
1921                              const register double tmp13_1 = tmp8_0*w1;                              const double tmp13_1 = tmp8_0*w1;
1922                              const register double tmp36_1 = tmp16_0*w1;                              const double tmp36_1 = tmp16_0*w1;
1923                              const register double tmp41_1 = A_01_3*w2;                              const double tmp41_1 = A_01_3*w2;
1924                              const register double tmp30_1 = tmp12_0*w13;                              const double tmp30_1 = tmp12_0*w13;
1925                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1926                              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(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
1927                              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(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
# Line 1659  void Rectangle::assemblePDESingle(Paso_S Line 1939  void Rectangle::assemblePDESingle(Paso_S
1939                              EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;                              EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
1940                              EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1941                          } else { // constant data                          } else { // constant data
1942                              const register double A_00 = A_p[INDEX2(0,0,2)];                              const double A_00 = A_p[INDEX2(0,0,2)];
1943                              const register double A_10 = A_p[INDEX2(1,0,2)];                              const double A_10 = A_p[INDEX2(1,0,2)];
1944                              const register double A_01 = A_p[INDEX2(0,1,2)];                              const double A_01 = A_p[INDEX2(0,1,2)];
1945                              const register double A_11 = A_p[INDEX2(1,1,2)];                              const double A_11 = A_p[INDEX2(1,1,2)];
1946                              const register double tmp0_0 = A_01 + A_10;                              const double tmp0_0 = A_01 + A_10;
1947                              const register double tmp0_1 = A_00*w18;                              const double tmp0_1 = A_00*w18;
1948                              const register double tmp1_1 = A_01*w19;                              const double tmp1_1 = A_01*w19;
1949                              const register double tmp2_1 = A_10*w20;                              const double tmp2_1 = A_10*w20;
1950                              const register double tmp3_1 = A_11*w21;                              const double tmp3_1 = A_11*w21;
1951                              const register double tmp4_1 = A_00*w22;                              const double tmp4_1 = A_00*w22;
1952                              const register double tmp5_1 = tmp0_0*w19;                              const double tmp5_1 = tmp0_0*w19;
1953                              const register double tmp6_1 = A_11*w23;                              const double tmp6_1 = A_11*w23;
1954                              const register double tmp7_1 = A_11*w25;                              const double tmp7_1 = A_11*w25;
1955                              const register double tmp8_1 = A_00*w24;                              const double tmp8_1 = A_00*w24;
1956                              const register double tmp9_1 = tmp0_0*w20;                              const double tmp9_1 = tmp0_0*w20;
1957                              const register double tmp10_1 = A_01*w20;                              const double tmp10_1 = A_01*w20;
1958                              const register double tmp11_1 = A_11*w27;                              const double tmp11_1 = A_11*w27;
1959                              const register double tmp12_1 = A_00*w26;                              const double tmp12_1 = A_00*w26;
1960                              const register double tmp13_1 = A_10*w19;                              const double tmp13_1 = A_10*w19;
1961                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1962                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1963                              EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
# Line 1703  void Rectangle::assemblePDESingle(Paso_S Line 1983  void Rectangle::assemblePDESingle(Paso_S
1983                          add_EM_S=true;                          add_EM_S=true;
1984                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1985                          if (B.actsExpanded()) {                          if (B.actsExpanded()) {
1986                              const register double B_0_0 = B_p[INDEX2(0,0,2)];                              const double B_0_0 = B_p[INDEX2(0,0,2)];
1987                              const register double B_1_0 = B_p[INDEX2(1,0,2)];                              const double B_1_0 = B_p[INDEX2(1,0,2)];
1988                              const register double B_0_1 = B_p[INDEX2(0,1,2)];                              const double B_0_1 = B_p[INDEX2(0,1,2)];
1989                              const register double B_1_1 = B_p[INDEX2(1,1,2)];                              const double B_1_1 = B_p[INDEX2(1,1,2)];
1990                              const register double B_0_2 = B_p[INDEX2(0,2,2)];                              const double B_0_2 = B_p[INDEX2(0,2,2)];
1991                              const register double B_1_2 = B_p[INDEX2(1,2,2)];                              const double B_1_2 = B_p[INDEX2(1,2,2)];
1992                              const register double B_0_3 = B_p[INDEX2(0,3,2)];                              const double B_0_3 = B_p[INDEX2(0,3,2)];
1993                              const register double B_1_3 = B_p[INDEX2(1,3,2)];                              const double B_1_3 = B_p[INDEX2(1,3,2)];
1994                              const register double tmp0_0 = B_1_0 + B_1_1;                              const double tmp0_0 = B_1_0 + B_1_1;
1995                              const register double tmp1_0 = B_1_2 + B_1_3;                              const double tmp1_0 = B_1_2 + B_1_3;
1996                              const register double tmp2_0 = B_0_1 + B_0_3;                              const double tmp2_0 = B_0_1 + B_0_3;
1997                              const register double tmp3_0 = B_0_0 + B_0_2;                              const double tmp3_0 = B_0_0 + B_0_2;
1998                              const register double tmp63_1 = B_1_1*w42;                              const double tmp63_1 = B_1_1*w42;
1999                              const register double tmp79_1 = B_1_1*w40;                              const double tmp79_1 = B_1_1*w40;
2000                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
2001                              const register double tmp8_1 = tmp0_0*w32;                              const double tmp8_1 = tmp0_0*w32;
2002                              const register double tmp71_1 = B_0_1*w34;                              const double tmp71_1 = B_0_1*w34;
2003                              const register double tmp19_1 = B_0_3*w31;                              const double tmp19_1 = B_0_3*w31;
2004                              const register double tmp15_1 = B_0_3*w34;                              const double tmp15_1 = B_0_3*w34;
2005                              const register double tmp9_1 = tmp3_0*w34;                              const double tmp9_1 = tmp3_0*w34;
2006                              const register double tmp35_1 = B_1_0*w36;                              const double tmp35_1 = B_1_0*w36;
2007                              const register double tmp66_1 = B_0_3*w28;                              const double tmp66_1 = B_0_3*w28;
2008                              const register double tmp28_1 = B_1_0*w42;                              const double tmp28_1 = B_1_0*w42;
2009                              const register double tmp22_1 = B_1_0*w40;                              const double tmp22_1 = B_1_0*w40;
2010                              const register double tmp16_1 = B_1_2*w29;                              const double tmp16_1 = B_1_2*w29;
2011                              const register double tmp6_1 = tmp2_0*w35;                              const double tmp6_1 = tmp2_0*w35;
2012                              const register double tmp55_1 = B_1_3*w40;                              const double tmp55_1 = B_1_3*w40;
2013                              const register double tmp50_1 = B_1_3*w42;                              const double tmp50_1 = B_1_3*w42;
2014                              const register double tmp7_1 = tmp1_0*w29;                              const double tmp7_1 = tmp1_0*w29;
2015                              const register double tmp1_1 = tmp1_0*w32;                              const double tmp1_1 = tmp1_0*w32;
2016                              const register double tmp57_1 = B_0_3*w30;                              const double tmp57_1 = B_0_3*w30;
2017                              const register double tmp18_1 = B_1_1*w32;                              const double tmp18_1 = B_1_1*w32;
2018                              const register double tmp53_1 = B_1_0*w41;                              const double tmp53_1 = B_1_0*w41;
2019                              const register double tmp61_1 = B_1_3*w36;                              const double tmp61_1 = B_1_3*w36;
2020                              const register double tmp27_1 = B_0_3*w38;                              const double tmp27_1 = B_0_3*w38;
2021                              const register double tmp64_1 = B_0_2*w30;                              const double tmp64_1 = B_0_2*w30;
2022                              const register double tmp76_1 = B_0_1*w38;                              const double tmp76_1 = B_0_1*w38;
2023                              const register double tmp39_1 = tmp2_0*w34;                              const double tmp39_1 = tmp2_0*w34;
2024                              const register double tmp62_1 = B_0_1*w31;                              const double tmp62_1 = B_0_1*w31;
2025                              const register double tmp56_1 = B_0_0*w31;                              const double tmp56_1 = B_0_0*w31;
2026                              const register double tmp49_1 = B_1_1*w36;                              const double tmp49_1 = B_1_1*w36;
2027                              const register double tmp2_1 = B_0_2*w31;                              const double tmp2_1 = B_0_2*w31;
2028                              const register double tmp23_1 = B_0_2*w33;                              const double tmp23_1 = B_0_2*w33;
2029                              const register double tmp38_1 = B_1_1*w43;                              const double tmp38_1 = B_1_1*w43;
2030                              const register double tmp74_1 = B_1_2*w41;                              const double tmp74_1 = B_1_2*w41;
2031                              const register double tmp43_1 = B_1_1*w41;                              const double tmp43_1 = B_1_1*w41;
2032                              const register double tmp58_1 = B_0_2*w28;                              const double tmp58_1 = B_0_2*w28;
2033                              const register double tmp67_1 = B_0_0*w33;                              const double tmp67_1 = B_0_0*w33;
2034                              const register double tmp33_1 = tmp0_0*w39;                              const double tmp33_1 = tmp0_0*w39;
2035                              const register double tmp4_1 = B_0_0*w28;                              const double tmp4_1 = B_0_0*w28;
2036                              const register double tmp20_1 = B_0_0*w30;                              const double tmp20_1 = B_0_0*w30;
2037                              const register double tmp13_1 = B_0_2*w38;                              const double tmp13_1 = B_0_2*w38;
2038                              const register double tmp65_1 = B_1_2*w43;                              const double tmp65_1 = B_1_2*w43;
2039                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
2040                              const register double tmp41_1 = tmp3_0*w33;                              const double tmp41_1 = tmp3_0*w33;
2041                              const register double tmp73_1 = B_0_2*w37;                              const double tmp73_1 = B_0_2*w37;
2042                              const register double tmp69_1 = B_0_0*w38;                              const double tmp69_1 = B_0_0*w38;
2043                              const register double tmp48_1 = B_1_2*w39;                              const double tmp48_1 = B_1_2*w39;
2044                              const register double tmp59_1 = B_0_1*w33;                              const double tmp59_1 = B_0_1*w33;
2045                              const register double tmp17_1 = B_1_3*w41;                              const double tmp17_1 = B_1_3*w41;
2046                              const register double tmp5_1 = B_0_3*w33;                              const double tmp5_1 = B_0_3*w33;
2047                              const register double tmp3_1 = B_0_1*w30;                              const double tmp3_1 = B_0_1*w30;
2048                              const register double tmp21_1 = B_0_1*w28;                              const double tmp21_1 = B_0_1*w28;
2049                              const register double tmp42_1 = B_1_0*w29;                              const double tmp42_1 = B_1_0*w29;
2050                              const register double tmp54_1 = B_1_2*w32;                              const double tmp54_1 = B_1_2*w32;
2051                              const register double tmp60_1 = B_1_0*w39;                              const double tmp60_1 = B_1_0*w39;
2052                              const register double tmp32_1 = tmp1_0*w36;                              const double tmp32_1 = tmp1_0*w36;
2053                              const register double tmp10_1 = B_0_1*w37;                              const double tmp10_1 = B_0_1*w37;
2054                              const register double tmp14_1 = B_0_0*w35;                              const double tmp14_1 = B_0_0*w35;
2055                              const register double tmp29_1 = B_0_1*w35;                              const double tmp29_1 = B_0_1*w35;
2056                              const register double tmp26_1 = B_1_2*w36;                              const double tmp26_1 = B_1_2*w36;
2057                              const register double tmp30_1 = B_1_3*w43;                              const double tmp30_1 = B_1_3*w43;
2058                              const register double tmp70_1 = B_0_2*w35;                              const double tmp70_1 = B_0_2*w35;
2059                              const register double tmp34_1 = B_1_3*w39;                              const double tmp34_1 = B_1_3*w39;
2060                              const register double tmp51_1 = B_1_0*w43;                              const double tmp51_1 = B_1_0*w43;
2061                              const register double tmp31_1 = B_0_2*w34;                              const double tmp31_1 = B_0_2*w34;
2062                              const register double tmp45_1 = tmp3_0*w28;                              const double tmp45_1 = tmp3_0*w28;
2063                              const register double tmp11_1 = tmp1_0*w39;                              const double tmp11_1 = tmp1_0*w39;
2064                              const register double tmp52_1 = B_1_1*w29;                              const double tmp52_1 = B_1_1*w29;
2065                              const register double tmp44_1 = B_1_3*w32;                              const double tmp44_1 = B_1_3*w32;
2066                              const register double tmp25_1 = B_1_1*w39;                              const double tmp25_1 = B_1_1*w39;
2067                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
2068                              const register double tmp72_1 = B_1_3*w29;                              const double tmp72_1 = B_1_3*w29;
2069                              const register double tmp40_1 = tmp2_0*w28;                              const double tmp40_1 = tmp2_0*w28;
2070                              const register double tmp46_1 = B_1_2*w40;                              const double tmp46_1 = B_1_2*w40;
2071                              const register double tmp36_1 = B_1_2*w42;                              const double tmp36_1 = B_1_2*w42;
2072                              const register double tmp24_1 = B_0_0*w37;                              const double tmp24_1 = B_0_0*w37;
2073                              const register double tmp77_1 = B_0_3*w35;                              const double tmp77_1 = B_0_3*w35;
2074                              const register double tmp68_1 = B_0_3*w37;                              const double tmp68_1 = B_0_3*w37;
2075                              const register double tmp78_1 = B_0_0*w34;                              const double tmp78_1 = B_0_0*w34;
2076                              const register double tmp12_1 = tmp0_0*w36;                              const double tmp12_1 = tmp0_0*w36;
2077                              const register double tmp75_1 = B_1_0*w32;                              const double tmp75_1 = B_1_0*w32;
2078                              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(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
2079                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2080                              EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
# Line 1812  void Rectangle::assemblePDESingle(Paso_S Line 2092  void Rectangle::assemblePDESingle(Paso_S
2092                              EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;                              EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2093                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              EM_S[INDEX2(3,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 register double B_0 = B_p[0];                              const double B_0 = B_p[0];
2096                              const register double B_1 = B_p[1];                              const double B_1 = B_p[1];
2097                              const register double tmp0_1 = B_0*w44;                              const double tmp0_1 = B_0*w44;
2098                              const register double tmp1_1 = B_1*w45;                              const double tmp1_1 = B_1*w45;
2099                              const register double tmp2_1 = B_0*w46;                              const double tmp2_1 = B_0*w46;
2100                              const register double tmp3_1 = B_0*w47;                              const double tmp3_1 = B_0*w47;
2101                              const register double tmp4_1 = B_1*w48;                              const double tmp4_1 = B_1*w48;
2102                              const register double tmp5_1 = B_1*w49;                              const double tmp5_1 = B_1*w49;
2103                              const register double tmp6_1 = B_1*w50;                              const double tmp6_1 = B_1*w50;
2104                              const register double tmp7_1 = B_0*w51;                              const double tmp7_1 = B_0*w51;
2105                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
2106                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2107                              EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;                              EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
# Line 1847  void Rectangle::assemblePDESingle(Paso_S Line 2127  void Rectangle::assemblePDESingle(Paso_S
2127                          add_EM_S=true;                          add_EM_S=true;
2128                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2129                          if (C.actsExpanded()) {                          if (C.actsExpanded()) {
2130                              const register double C_0_0 = C_p[INDEX2(0,0,2)];                              const double C_0_0 = C_p[INDEX2(0,0,2)];
2131                              const register double C_1_0 = C_p[INDEX2(1,0,2)];                              const double C_1_0 = C_p[INDEX2(1,0,2)];
2132                              const register double C_0_1 = C_p[INDEX2(0,1,2)];                              const double C_0_1 = C_p[INDEX2(0,1,2)];
2133                              const register double C_1_1 = C_p[INDEX2(1,1,2)];                              const double C_1_1 = C_p[INDEX2(1,1,2)];
2134                              const register double C_0_2 = C_p[INDEX2(0,2,2)];                              const double C_0_2 = C_p[INDEX2(0,2,2)];
2135                              const register double C_1_2 = C_p[INDEX2(1,2,2)];                              const double C_1_2 = C_p[INDEX2(1,2,2)];
2136                              const register double C_0_3 = C_p[INDEX2(0,3,2)];                              const double C_0_3 = C_p[INDEX2(0,3,2)];
2137                              const register double C_1_3 = C_p[INDEX2(1,3,2)];                              const double C_1_3 = C_p[INDEX2(1,3,2)];
2138                              const register double tmp0_0 = C_1_0 + C_1_1;                              const double tmp0_0 = C_1_0 + C_1_1;
2139                              const register double tmp1_0 = C_1_2 + C_1_3;                              const double tmp1_0 = C_1_2 + C_1_3;
2140                              const register double tmp2_0 = C_0_1 + C_0_3;                              const double tmp2_0 = C_0_1 + C_0_3;
2141                              const register double tmp3_0 = C_0_0 + C_0_2;                              const double tmp3_0 = C_0_0 + C_0_2;
2142                              const register double tmp64_1 = C_0_2*w30;                              const double tmp64_1 = C_0_2*w30;
2143                              const register double tmp14_1 = C_0_2*w28;                              const double tmp14_1 = C_0_2*w28;
2144                              const register double tmp19_1 = C_0_3*w31;                              const double tmp19_1 = C_0_3*w31;
2145                              const register double tmp22_1 = C_1_0*w40;                              const double tmp22_1 = C_1_0*w40;
2146                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
2147                              const register double tmp29_1 = C_0_1*w35;                              const double tmp29_1 = C_0_1*w35;
2148                              const register double tmp73_1 = C_0_2*w37;                              const double tmp73_1 = C_0_2*w37;
2149                              const register double tmp74_1 = C_1_2*w41;                              const double tmp74_1 = C_1_2*w41;
2150                              const register double tmp52_1 = C_1_3*w39;                              const double tmp52_1 = C_1_3*w39;
2151                              const register double tmp25_1 = C_1_1*w39;                              const double tmp25_1 = C_1_1*w39;
2152                              const register double tmp62_1 = C_0_1*w31;                              const double tmp62_1 = C_0_1*w31;
2153                              const register double tmp79_1 = C_1_1*w40;                              const double tmp79_1 = C_1_1*w40;
2154                              const register double tmp43_1 = C_1_1*w36;                              const double tmp43_1 = C_1_1*w36;
2155                              const register double tmp27_1 = C_0_3*w38;                              const double tmp27_1 = C_0_3*w38;
2156                              const register double tmp28_1 = C_1_0*w42;                              const double tmp28_1 = C_1_0*w42;
2157                              const register double tmp63_1 = C_1_1*w42;                              const double tmp63_1 = C_1_1*w42;
2158                              const register double tmp59_1 = C_0_3*w34;                              const double tmp59_1 = C_0_3*w34;
2159                              const register double tmp72_1 = C_1_3*w29;                              const double tmp72_1 = C_1_3*w29;
2160                              const register double tmp40_1 = tmp2_0*w35;                              const double tmp40_1 = tmp2_0*w35;
2161                              const register double tmp13_1 = C_0_3*w30;                              const double tmp13_1 = C_0_3*w30;
2162                              const register double tmp51_1 = C_1_2*w40;                              const double tmp51_1 = C_1_2*w40;
2163                              const register double tmp54_1 = C_1_2*w42;                              const double tmp54_1 = C_1_2*w42;
2164                              const register double tmp12_1 = C_0_0*w31;                              const double tmp12_1 = C_0_0*w31;
2165                              const register double tmp2_1 = tmp1_0*w32;                              const double tmp2_1 = tmp1_0*w32;
2166                              const register double tmp68_1 = C_0_2*w31;                              const double tmp68_1 = C_0_2*w31;
2167                              const register double tmp75_1 = C_1_0*w32;                              const double tmp75_1 = C_1_0*w32;
2168                              const register double tmp49_1 = C_1_1*w41;                              const double tmp49_1 = C_1_1*w41;
2169                              const register double tmp4_1 = C_0_2*w35;                              const double tmp4_1 = C_0_2*w35;
2170                              const register double tmp66_1 = C_0_3*w28;                              const double tmp66_1 = C_0_3*w28;
2171                              const register double tmp56_1 = C_0_1*w37;                              const double tmp56_1 = C_0_1*w37;
2172                              const register double tmp5_1 = C_0_1*w34;                              const double tmp5_1 = C_0_1*w34;
2173                              const register double tmp38_1 = tmp2_0*w34;                              const double tmp38_1 = tmp2_0*w34;
2174                              const register double tmp76_1 = C_0_1*w38;                              const double tmp76_1 = C_0_1*w38;
2175                              const register double tmp21_1 = C_0_1*w28;                              const double tmp21_1 = C_0_1*w28;
2176                              const register double tmp69_1 = C_0_1*w30;                              const double tmp69_1 = C_0_1*w30;
2177                              const register double tmp53_1 = C_1_0*w36;                              const double tmp53_1 = C_1_0*w36;
2178                              const register double tmp42_1 = C_1_2*w39;                              const double tmp42_1 = C_1_2*w39;
2179                              const register double tmp32_1 = tmp1_0*w29;                              const double tmp32_1 = tmp1_0*w29;
2180                              const register double tmp45_1 = C_1_0*w43;                              const double tmp45_1 = C_1_0*w43;
2181                              const register double tmp33_1 = tmp0_0*w32;                              const double tmp33_1 = tmp0_0*w32;
2182                              const register double tmp35_1 = C_1_0*w41;                              const double tmp35_1 = C_1_0*w41;
2183                              const register double tmp26_1 = C_1_2*w36;                              const double tmp26_1 = C_1_2*w36;
2184                              const register double tmp67_1 = C_0_0*w33;                              const double tmp67_1 = C_0_0*w33;
2185                              const register double tmp31_1 = C_0_2*w34;                              const double tmp31_1 = C_0_2*w34;
2186                              const register double tmp20_1 = C_0_0*w30;                              const double tmp20_1 = C_0_0*w30;
2187                              const register double tmp70_1 = C_0_0*w28;                              const double tmp70_1 = C_0_0*w28;
2188                              const register double tmp8_1 = tmp0_0*w39;                              const double tmp8_1 = tmp0_0*w39;
2189                              const register double tmp30_1 = C_1_3*w43;                              const double tmp30_1 = C_1_3*w43;
2190                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
2191                              const register double tmp17_1 = C_1_3*w41;                              const double tmp17_1 = C_1_3*w41;
2192                              const register double tmp58_1 = C_0_0*w35;                              const double tmp58_1 = C_0_0*w35;
2193                              const register double tmp9_1 = tmp3_0*w33;                              const double tmp9_1 = tmp3_0*w33;
2194                              const register double tmp61_1 = C_1_3*w36;                              const double tmp61_1 = C_1_3*w36;
2195                              const register double tmp41_1 = tmp3_0*w34;                              const double tmp41_1 = tmp3_0*w34;
2196                              const register double tmp50_1 = C_1_3*w32;                              const double tmp50_1 = C_1_3*w32;
2197                              const register double tmp18_1 = C_1_1*w32;                              const double tmp18_1 = C_1_1*w32;
2198                              const register double tmp6_1 = tmp1_0*w36;                              const double tmp6_1 = tmp1_0*w36;
2199                              const register double tmp3_1 = C_0_0*w38;                              const double tmp3_1 = C_0_0*w38;
2200                              const register double tmp34_1 = C_1_1*w29;                              const double tmp34_1 = C_1_1*w29;
2201                              const register double tmp77_1 = C_0_3*w35;                              const double tmp77_1 = C_0_3*w35;
2202                              const register double tmp65_1 = C_1_2*w43;                              const double tmp65_1 = C_1_2*w43;
2203                              const register double tmp71_1 = C_0_3*w33;                              const double tmp71_1 = C_0_3*w33;
2204                              const register double tmp55_1 = C_1_1*w43;                              const double tmp55_1 = C_1_1*w43;
2205                              const register double tmp46_1 = tmp3_0*w28;                              const double tmp46_1 = tmp3_0*w28;
2206                              const register double tmp24_1 = C_0_0*w37;                              const double tmp24_1 = C_0_0*w37;
2207                              const register double tmp10_1 = tmp1_0*w39;                              const double tmp10_1 = tmp1_0*w39;
2208                              const register double tmp48_1 = C_1_0*w29;                              const double tmp48_1 = C_1_0*w29;
2209                              const register double tmp15_1 = C_0_1*w33;                              const double tmp15_1 = C_0_1*w33;
2210                              const register double tmp36_1 = C_1_2*w32;                              const double tmp36_1 = C_1_2*w32;
2211                              const register double tmp60_1 = C_1_0*w39;                              const double tmp60_1 = C_1_0*w39;
2212                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
2213                              const register double tmp16_1 = C_1_2*w29;                              const double tmp16_1 = C_1_2*w29;
2214                              const register double tmp1_1 = C_0_3*w37;                              const double tmp1_1 = C_0_3*w37;
2215                              const register double tmp7_1 = tmp2_0*w28;                              const double tmp7_1 = tmp2_0*w28;
2216                              const register double tmp39_1 = C_1_3*w40;                              const double tmp39_1 = C_1_3*w40;
2217                              const register double tmp44_1 = C_1_3*w42;                              const double tmp44_1 = C_1_3*w42;
2218                              const register double tmp57_1 = C_0_2*w38;                              const double tmp57_1 = C_0_2*w38;
2219                              const register double tmp78_1 = C_0_0*w34;                              const double tmp78_1 = C_0_0*w34;
2220                              const register double tmp11_1 = tmp0_0*w36;                              const double tmp11_1 = tmp0_0*w36;
2221                              const register double tmp23_1 = C_0_2*w33;                              const double tmp23_1 = C_0_2*w33;
2222                              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(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
2223                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2224                              EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
# Line 1956  void Rectangle::assemblePDESingle(Paso_S Line 2236  void Rectangle::assemblePDESingle(Paso_S
2236                              EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;                              EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2237                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2238                          } else { // constant data                          } else { // constant data
2239                              const register double C_0 = C_p[0];                              const double C_0 = C_p[0];
2240                              const register double C_1 = C_p[1];                              const double C_1 = C_p[1];
2241                              const register double tmp0_1 = C_0*w47;                              const double tmp0_1 = C_0*w47;
2242                              const register double tmp1_1 = C_1*w45;                              const double tmp1_1 = C_1*w45;
2243                              const register double tmp2_1 = C_1*w48;                              const double tmp2_1 = C_1*w48;
2244                              const register double tmp3_1 = C_0*w51;                              const double tmp3_1 = C_0*w51;
2245                              const register double tmp4_1 = C_0*w44;                              const double tmp4_1 = C_0*w44;
2246                              const register double tmp5_1 = C_1*w49;                              const double tmp5_1 = C_1*w49;
2247                              const register double tmp6_1 = C_1*w50;                              const double tmp6_1 = C_1*w50;
2248                              const register double tmp7_1 = C_0*w46;                              const double tmp7_1 = C_0*w46;
2249                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
2250                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
2251                              EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
# Line 1991  void Rectangle::assemblePDESingle(Paso_S Line 2271  void Rectangle::assemblePDESingle(Paso_S
2271                          add_EM_S=true;                          add_EM_S=true;
2272                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2273                          if (D.actsExpanded()) {                          if (D.actsExpanded()) {
2274                              const register double D_0 = D_p[0];                              const double D_0 = D_p[0];
2275                              const register double D_1 = D_p[1];                              const double D_1 = D_p[1];
2276                              const register double D_2 = D_p[2];                              const double D_2 = D_p[2];
2277                              const register double D_3 = D_p[3];                              const double D_3 = D_p[3];
2278                              const register double tmp0_0 = D_0 + D_1;                              const double tmp0_0 = D_0 + D_1;
2279                              const register double tmp1_0 = D_2 + D_3;                              const double tmp1_0 = D_2 + D_3;
2280                              const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;                              const double tmp2_0 = D_0 + D_1 + D_2 + D_3;
2281                              const register double tmp3_0 = D_1 + D_2;                              const double tmp3_0 = D_1 + D_2;
2282                              const register double tmp4_0 = D_1 + D_3;                              const double tmp4_0 = D_1 + D_3;
2283                              const register double tmp5_0 = D_0 + D_2;                              const double tmp5_0 = D_0 + D_2;
2284                              const register double tmp6_0 = D_0 + D_3;                              const double tmp6_0 = D_0 + D_3;
2285                              const register double tmp0_1 = tmp0_0*w52;                              const double tmp0_1 = tmp0_0*w52;
2286                              const register double tmp1_1 = tmp1_0*w53;                              const double tmp1_1 = tmp1_0*w53;
2287                              const register double tmp2_1 = tmp2_0*w54;                              const double tmp2_1 = tmp2_0*w54;
2288                              const register double tmp3_1 = tmp1_0*w52;                              const double tmp3_1 = tmp1_0*w52;
2289                              const register double tmp4_1 = tmp0_0*w53;                              const double tmp4_1 = tmp0_0*w53;
2290                              const register double tmp5_1 = tmp3_0*w54;                              const double tmp5_1 = tmp3_0*w54;
2291                              const register double tmp6_1 = D_0*w55;                              const double tmp6_1 = D_0*w55;
2292                              const register double tmp7_1 = D_3*w56;                              const double tmp7_1 = D_3*w56;
2293                              const register double tmp8_1 = D_3*w55;                              const double tmp8_1 = D_3*w55;
2294                              const register double tmp9_1 = D_0*w56;                              const double tmp9_1 = D_0*w56;
2295                              const register double tmp10_1 = tmp4_0*w52;                              const double tmp10_1 = tmp4_0*w52;
2296                              const register double tmp11_1 = tmp5_0*w53;                              const double tmp11_1 = tmp5_0*w53;
2297                              const register double tmp12_1 = tmp5_0*w52;                              const double tmp12_1 = tmp5_0*w52;
2298                              const register double tmp13_1 = tmp4_0*w53;                              const double tmp13_1 = tmp4_0*w53;
2299                              const register double tmp14_1 = tmp6_0*w54;                              const double tmp14_1 = tmp6_0*w54;
2300                              const register double tmp15_1 = D_2*w55;                              const double tmp15_1 = D_2*w55;
2301                              const register double tmp16_1 = D_1*w56;                              const double tmp16_1 = D_1*w56;
2302                              const register double tmp17_1 = D_1*w55;                              const double tmp17_1 = D_1*w55;
2303                              const register double tmp18_1 = D_2*w56;                              const double tmp18_1 = D_2*w56;
2304                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2305                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2306                              EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;                              EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
# Line 2038  void Rectangle::assemblePDESingle(Paso_S Line 2318  void Rectangle::assemblePDESingle(Paso_S
2318                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2319                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2320                          } else { // constant data                          } else { // constant data
2321                              const register double tmp0_1 = D_p[0]*w57;                              const double tmp0_1 = D_p[0]*w57;
2322                              const register double tmp1_1 = D_p[0]*w58;                              const double tmp1_1 = D_p[0]*w58;
2323                              const register double tmp2_1 = D_p[0]*w59;                              const double tmp2_1 = D_p[0]*w59;
2324                              EM_S[INDEX2(0,0,4)]+=tmp2_1;                              EM_S[INDEX2(0,0,4)]+=tmp2_1;
2325                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
2326                              EM_S[INDEX2(2,0,4)]+=tmp0_1;                              EM_S[INDEX2(2,0,4)]+=tmp0_1;
# Line 2066  void Rectangle::assemblePDESingle(Paso_S Line 2346  void Rectangle::assemblePDESingle(Paso_S
2346                          add_EM_F=true;                          add_EM_F=true;
2347                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2348                          if (X.actsExpanded()) {                          if (X.actsExpanded()) {
2349                              const register double X_0_0 = X_p[INDEX2(0,0,2)];                              const double X_0_0 = X_p[INDEX2(0,0,2)];
2350                              const register double X_1_0 = X_p[INDEX2(1,0,2)];                              const double X_1_0 = X_p[INDEX2(1,0,2)];
2351                              const register double X_0_1 = X_p[INDEX2(0,1,2)];                              const double X_0_1 = X_p[INDEX2(0,1,2)];
2352                              const register double X_1_1 = X_p[INDEX2(1,1,2)];                              const double X_1_1 = X_p[INDEX2(1,1,2)];
2353                              const register double X_0_2 = X_p[INDEX2(0,2,2)];                              const double X_0_2 = X_p[INDEX2(0,2,2)];
2354                              const register double X_1_2 = X_p[INDEX2(1,2,2)];                              const double X_1_2 = X_p[INDEX2(1,2,2)];
2355                              const register double X_0_3 = X_p[INDEX2(0,3,2)];                              const double X_0_3 = X_p[INDEX2(0,3,2)];
2356                              const register double X_1_3 = X_p[INDEX2(1,3,2)];                              const double X_1_3 = X_p[INDEX2(1,3,2)];
2357                              const register double tmp0_0 = X_0_2 + X_0_3;                              const double tmp0_0 = X_0_2 + X_0_3;
2358                              const register double tmp1_0 = X_1_1 + X_1_3;                              const double tmp1_0 = X_1_1 + X_1_3;
2359                              const register double tmp2_0 = X_1_0 + X_1_2;                              const double tmp2_0 = X_1_0 + X_1_2;
2360                              const register double tmp3_0 = X_0_0 + X_0_1;                              const double tmp3_0 = X_0_0 + X_0_1;
2361                              const register double tmp0_1 = tmp0_0*w63;                              const double tmp0_1 = tmp0_0*w63;
2362                              const register double tmp1_1 = tmp1_0*w62;                              const double tmp1_1 = tmp1_0*w62;
2363                              const register double tmp2_1 = tmp2_0*w61;                              const double tmp2_1 = tmp2_0*w61;
2364                              const register double tmp3_1 = tmp3_0*w60;                              const double tmp3_1 = tmp3_0*w60;
2365                              const register double tmp4_1 = tmp0_0*w65;                              const double tmp4_1 = tmp0_0*w65;
2366                              const register double tmp5_1 = tmp3_0*w64;                              const double tmp5_1 = tmp3_0*w64;
2367                              const register double tmp6_1 = tmp2_0*w62;                              const double tmp6_1 = tmp2_0*w62;
2368                              const register double tmp7_1 = tmp1_0*w61;                              const double tmp7_1 = tmp1_0*w61;
2369                              const register double tmp8_1 = tmp2_0*w66;                              const double tmp8_1 = tmp2_0*w66;
2370                              const register double tmp9_1 = tmp3_0*w63;                              const double tmp9_1 = tmp3_0*w63;
2371                              const register double tmp10_1 = tmp0_0*w60;                              const double tmp10_1 = tmp0_0*w60;
2372                              const register double tmp11_1 = tmp1_0*w67;                              const double tmp11_1 = tmp1_0*w67;
2373                              const register double tmp12_1 = tmp1_0*w66;                              const double tmp12_1 = tmp1_0*w66;
2374                              const register double tmp13_1 = tmp3_0*w65;                              const double tmp13_1 = tmp3_0*w65;
2375                              const register double tmp14_1 = tmp0_0*w64;                              const double tmp14_1 = tmp0_0*w64;
2376                              const register double tmp15_1 = tmp2_0*w67;                              const double tmp15_1 = tmp2_0*w67;
2377                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2378                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2379                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2380                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2381                          } else { // constant data                          } else { // constant data
2382                              const register double tmp0_1 = X_p[1]*w69;                              const double tmp0_1 = X_p[1]*w69;
2383                              const register double tmp1_1 = X_p[0]*w68;                              const double tmp1_1 = X_p[0]*w68;
2384                              const register double tmp2_1 = X_p[0]*w70;                              const double tmp2_1 = X_p[0]*w70;
2385                              const register double tmp3_1 = X_p[1]*w71;                              const double tmp3_1 = X_p[1]*w71;
2386                              EM_F[0]+=tmp0_1 + tmp1_1;                              EM_F[0]+=tmp0_1 + tmp1_1;
2387                              EM_F[1]+=tmp0_1 + tmp2_1;                              EM_F[1]+=tmp0_1 + tmp2_1;
2388                              EM_F[2]+=tmp1_1 + tmp3_1;                              EM_F[2]+=tmp1_1 + tmp3_1;
# Line 2116  void Rectangle::assemblePDESingle(Paso_S Line 2396  void Rectangle::assemblePDESingle(Paso_S
2396                          add_EM_F=true;                          add_EM_F=true;
2397                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2398                          if (Y.actsExpanded()) {                          if (Y.actsExpanded()) {
2399                              const register double Y_0 = Y_p[0];                              const double Y_0 = Y_p[0];
2400                              const register double Y_1 = Y_p[1];                              const double Y_1 = Y_p[1];
2401                              const register double Y_2 = Y_p[2];                              const double Y_2 = Y_p[2];
2402                              const register double Y_3 = Y_p[3];                              const double Y_3 = Y_p[3];
2403                              const register double tmp0_0 = Y_1 + Y_2;                              const double tmp0_0 = Y_1 + Y_2;
2404                              const register double tmp1_0 = Y_0 + Y_3;                              const double tmp1_0 = Y_0 + Y_3;
2405                              const register double tmp0_1 = Y_0*w72;                              const double tmp0_1 = Y_0*w72;
2406                              const register double tmp1_1 = tmp0_0*w73;                              const double tmp1_1 = tmp0_0*w73;
2407                              const register double tmp2_1 = Y_3*w74;                              const double tmp2_1 = Y_3*w74;
2408                              const register double tmp3_1 = Y_1*w72;                              const double tmp3_1 = Y_1*w72;
2409                              const register double tmp4_1 = tmp1_0*w73;                              const double tmp4_1 = tmp1_0*w73;
2410                              const register double tmp5_1 = Y_2*w74;                              const double tmp5_1 = Y_2*w74;
2411                              const register double tmp6_1 = Y_2*w72;                              const double tmp6_1 = Y_2*w72;
2412                              const register double tmp7_1 = Y_1*w74;                              const double tmp7_1 = Y_1*w74;
2413                              const register double tmp8_1 = Y_3*w72;                              const double tmp8_1 = Y_3*w72;
2414                              const register double tmp9_1 = Y_0*w74;                              const double tmp9_1 = Y_0*w74;
2415                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2416                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2417                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2418                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2419                          } else { // constant data                          } else { // constant data
2420                              const register double tmp0_1 = Y_p[0]*w75;                              const double tmp0_1 = Y_p[0]*w75;
2421                              EM_F[0]+=tmp0_1;                              EM_F[0]+=tmp0_1;
2422                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
2423                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
# Line 2147  void Rectangle::assemblePDESingle(Paso_S Line 2427  void Rectangle::assemblePDESingle(Paso_S
2427    
2428                      // 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)
2429                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
2430                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
                     rowIndex.push_back(m_dofMap[firstNode]);  
                     rowIndex.push_back(m_dofMap[firstNode+1]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                     if (add_EM_F) {  
                         double *F_p=rhs.getSampleDataRW(0);  
                         for (index_t i=0; i<rowIndex.size(); i++) {  
                             if (rowIndex[i]<getNumDOF()) {  
                                 F_p[rowIndex[i]]+=EM_F[i];  
                             }  
                         }  
                     }  
                     if (add_EM_S) {  
                         addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);  
                     }  
   
2431                  } // end k0 loop                  } // end k0 loop
2432              } // end k1 loop              } // end k1 loop
2433          } // end of colouring          } // end of colouring
# Line 2174  void Rectangle::assemblePDESingle(Paso_S Line 2438  void Rectangle::assemblePDESingle(Paso_S
2438  void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2439          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2440          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2441          const escript::Data& X, const escript::Data& Y,          const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
2442  {  {
2443      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
2444      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 2214  void Rectangle::assemblePDESingleReduced Line 2477  void Rectangle::assemblePDESingleReduced
2477                      if (!A.isEmpty()) {                      if (!A.isEmpty()) {
2478                          add_EM_S=true;                          add_EM_S=true;
2479                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2480                          const register double A_00 = A_p[INDEX2(0,0,2)];                          const double A_00 = A_p[INDEX2(0,0,2)];
2481                          const register double A_10 = A_p[INDEX2(1,0,2)];                          const double A_10 = A_p[INDEX2(1,0,2)];
2482                          const register double A_01 = A_p[INDEX2(0,1,2)];                          const double A_01 = A_p[INDEX2(0,1,2)];
2483                          const register double A_11 = A_p[INDEX2(1,1,2)];                          const double A_11 = A_p[INDEX2(1,1,2)];
2484                          const register double tmp0_0 = A_01 + A_10;                          const double tmp0_0 = A_01 + A_10;
2485                          const register double tmp0_1 = A_11*w3;                          const double tmp0_1 = A_11*w3;
2486                          const register double tmp1_1 = A_00*w0;                          const double tmp1_1 = A_00*w0;
2487                          const register double tmp2_1 = A_01*w1;                          const double tmp2_1 = A_01*w1;
2488                          const register double tmp3_1 = A_10*w2;                          const double tmp3_1 = A_10*w2;
2489                          const register double tmp4_1 = tmp0_0*w1;                          const double tmp4_1 = tmp0_0*w1;
2490                          const register double tmp5_1 = A_11*w4;                          const double tmp5_1 = A_11*w4;
2491                          const register double tmp6_1 = A_00*w5;                          const double tmp6_1 = A_00*w5;
2492                          const register double tmp7_1 = tmp0_0*w2;                          const double tmp7_1 = tmp0_0*w2;
2493                          const register double tmp8_1 = A_10*w1;                          const double tmp8_1 = A_10*w1;
2494                          const register double tmp9_1 = A_01*w2;                          const double tmp9_1 = A_01*w2;
2495                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2496                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2497                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
# Line 2252  void Rectangle::assemblePDESingleReduced Line 2515  void Rectangle::assemblePDESingleReduced
2515                      if (!B.isEmpty()) {                      if (!B.isEmpty()) {
2516                          add_EM_S=true;                          add_EM_S=true;
2517                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2518                          const register double tmp2_1 = B_p[0]*w8;                          const double tmp2_1 = B_p[0]*w8;
2519                          const register double tmp0_1 = B_p[0]*w6;                          const double tmp0_1 = B_p[0]*w6;
2520                          const register double tmp3_1 = B_p[1]*w9;                          const double tmp3_1 = B_p[1]*w9;
2521                          const register double tmp1_1 = B_p[1]*w7;                          const double tmp1_1 = B_p[1]*w7;
2522                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2523                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2524                          EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
# Line 2279  void Rectangle::assemblePDESingleReduced Line 2542  void Rectangle::assemblePDESingleReduced
2542                      if (!C.isEmpty()) {                      if (!C.isEmpty()) {
2543                          add_EM_S=true;                          add_EM_S=true;
2544                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2545                          const register double tmp1_1 = C_p[1]*w7;                          const double tmp1_1 = C_p[1]*w7;
2546                          const register double tmp0_1 = C_p[0]*w8;                          const double tmp0_1 = C_p[0]*w8;
2547                          const register double tmp3_1 = C_p[0]*w6;                          const double tmp3_1 = C_p[0]*w6;
2548                          const register double tmp2_1 = C_p[1]*w9;                          const double tmp2_1 = C_p[1]*w9;
2549                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2550                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2551                          EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
# Line 2306  void Rectangle::assemblePDESingleReduced Line 2569  void Rectangle::assemblePDESingleReduced
2569                      if (!D.isEmpty()) {                      if (!D.isEmpty()) {
2570                          add_EM_S=true;                          add_EM_S=true;
2571                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2572                          const register double tmp0_1 = D_p[0]*w10;                          const double tmp0_1 = D_p[0]*w10;
2573                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1;
2574                          EM_S[INDEX2(1,0,4)]+=tmp0_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1;
2575                          EM_S[INDEX2(2,0,4)]+=tmp0_1;                          EM_S[INDEX2(2,0,4)]+=tmp0_1;
# Line 2330  void Rectangle::assemblePDESingleReduced Line 2593  void Rectangle::assemblePDESingleReduced
2593                      if (!X.isEmpty()) {                      if (!X.isEmpty()) {
2594                          add_EM_F=true;                          add_EM_F=true;
2595                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2596                          const register double tmp0_1 = X_p[0]*w11;                          const double tmp0_1 = X_p[0]*w11;
2597                          const register double tmp2_1 = X_p[0]*w13;                          const double tmp2_1 = X_p[0]*w13;
2598                          const register double tmp1_1 = X_p[1]*w12;                          const double tmp1_1 = X_p[1]*w12;
2599                          const register double tmp3_1 = X_p[1]*w14;                          const double tmp3_1 = X_p[1]*w14;
2600                          EM_F[0]+=tmp0_1 + tmp1_1;                          EM_F[0]+=tmp0_1 + tmp1_1;
2601                          EM_F[1]+=tmp1_1 + tmp2_1;                          EM_F[1]+=tmp1_1 + tmp2_1;
2602                          EM_F[2]+=tmp0_1 + tmp3_1;                          EM_F[2]+=tmp0_1 + tmp3_1;
# Line 2345  void Rectangle::assemblePDESingleReduced Line 2608  void Rectangle::assemblePDESingleReduced
2608                      if (!Y.isEmpty()) {                      if (!Y.isEmpty()) {
2609                          add_EM_F=true;                          add_EM_F=true;
2610                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2611                          const register double tmp0_1 = Y_p[0]*w15;                          const double tmp0_1 = Y_p[0]*w15;
2612                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0_1;
2613                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0_1;
2614                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
# Line 2354  void Rectangle::assemblePDESingleReduced Line 2617  void Rectangle::assemblePDESingleReduced
2617    
2618                      // 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)
2619                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
2620                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
                     rowIndex.push_back(m_dofMap[firstNode]);  
                     rowIndex.push_back(m_dofMap[firstNode+1]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                     if (add_EM_F) {  
                         double *F_p=rhs.getSampleDataRW(0);  
                         for (index_t i=0; i<rowIndex.size(); i++) {  
                             if (rowIndex[i]<getNumDOF()) {  
                                 F_p[rowIndex[i]]+=EM_F[i];  
                             }  
                         }  
                     }  
                     if (add_EM_S) {  
                         addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);  
                     }  
   
2621                  } // end k0 loop                  } // end k0 loop
2622              } // end k1 loop              } // end k1 loop
2623          } // end of colouring          } // end of colouring
# Line 2381  void Rectangle::assemblePDESingleReduced Line 2628  void Rectangle::assemblePDESingleReduced
2628  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2629          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2630          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2631          const escript::Data& X, const escript::Data& Y,          const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
2632  {  {
2633      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
2634      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 2492  void Rectangle::assemblePDESystem(Paso_S Line 2738  void Rectangle::assemblePDESystem(Paso_S
2738                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
2739                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2740                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2741                                      const register 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)];
2742                                      const register 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)];
2743                                      const register 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)];
2744                                      const register 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)];
2745                                      const register 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)];
2746                                      const register 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)];
2747                                      const register 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)];
2748                                      const register 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)];
2749                                      const register 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)];
2750                                      const register 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)];
2751                                      const register 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)];
2752                                      const register 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)];
2753                                      const register 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)];
2754                                      const register 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)];
2755                                      const register 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)];
2756                                      const register 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)];
2757                                      const register double tmp0_0 = A_01_0 + A_01_3;                                      const double tmp0_0 = A_01_0 + A_01_3;
2758                                      const register double tmp1_0 = A_00_0 + A_00_1;                                      const double tmp1_0 = A_00_0 + A_00_1;
2759                                      const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                                      const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
2760                                      const register double tmp3_0 = A_00_2 + A_00_3;                                      const double tmp3_0 = A_00_2 + A_00_3;
2761                                      const register double tmp4_0 = A_10_1 + A_10_2;                                      const double tmp4_0 = A_10_1 + A_10_2;
2762                                      const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                                      const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
2763                                      const register double tmp6_0 = A_01_3 + A_10_0;                                      const double tmp6_0 = A_01_3 + A_10_0;
2764                                      const register double tmp7_0 = A_01_0 + A_10_3;                                      const double tmp7_0 = A_01_0 + A_10_3;
2765                                      const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                                      const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
2766                                      const register double tmp9_0 = A_01_0 + A_10_0;                                      const double tmp9_0 = A_01_0 + A_10_0;
2767                                      const register double tmp10_0 = A_01_3 + A_10_3;                                      const double tmp10_0 = A_01_3 + A_10_3;
2768                                      const register double tmp11_0 = A_11_1 + A_11_3;                                      const double tmp11_0 = A_11_1 + A_11_3;
2769                                      const register double tmp12_0 = A_11_0 + A_11_2;                                      const double tmp12_0 = A_11_0 + A_11_2;
2770                                      const register double tmp13_0 = A_01_2 + A_10_1;                                      const double tmp13_0 = A_01_2 + A_10_1;
2771                                      const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                                      const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
2772                                      const register double tmp15_0 = A_01_1 + A_10_2;                                      const double tmp15_0 = A_01_1 + A_10_2;
2773                                      const register double tmp16_0 = A_10_0 + A_10_3;                                      const double tmp16_0 = A_10_0 + A_10_3;
2774                                      const register double tmp17_0 = A_01_1 + A_01_2;                                      const double tmp17_0 = A_01_1 + A_01_2;
2775                                      const register double tmp18_0 = A_01_1 + A_10_1;                                      const double tmp18_0 = A_01_1 + A_10_1;
2776                                      const register double tmp19_0 = A_01_2 + A_10_2;                                      const double tmp19_0 = A_01_2 + A_10_2;
2777                                      const register double tmp0_1 = A_10_3*w8;                                      const double tmp0_1 = A_10_3*w8;
2778                                      const register double tmp1_1 = tmp0_0*w1;                                      const double tmp1_1 = tmp0_0*w1;
2779                                      const register double tmp2_1 = A_01_1*w4;                                      const double tmp2_1 = A_01_1*w4;
2780                                      const register double tmp3_1 = tmp1_0*w0;                                      const double tmp3_1 = tmp1_0*w0;
2781                                      const register double tmp4_1 = A_01_2*w7;                                      const double tmp4_1 = A_01_2*w7;
2782                                      const register double tmp5_1 = tmp2_0*w3;                                      const double tmp5_1 = tmp2_0*w3;
2783                                      const register double tmp6_1 = tmp3_0*w6;                                      const double tmp6_1 = tmp3_0*w6;
2784                                      const register double tmp7_1 = A_10_0*w2;                                      const double tmp7_1 = A_10_0*w2;
2785                                      const register double tmp8_1 = tmp4_0*w5;                                      const double tmp8_1 = tmp4_0*w5;
2786                                      const register double tmp9_1 = tmp2_0*w10;                                      const double tmp9_1 = tmp2_0*w10;
2787                                      const register double tmp10_1 = tmp5_0*w9;                                      const double tmp10_1 = tmp5_0*w9;
2788                                      const register double tmp11_1 = tmp6_0*w7;                                      const double tmp11_1 = tmp6_0*w7;
2789                                      const register double tmp12_1 = tmp7_0*w4;                                      const double tmp12_1 = tmp7_0*w4;
2790                                      const register double tmp13_1 = tmp8_0*w1;                                      const double tmp13_1 = tmp8_0*w1;
2791                                      const register double tmp14_1 = A_10_0*w8;                                      const double tmp14_1 = A_10_0*w8;
2792                                      const register double tmp15_1 = A_01_2*w4;                                      const double tmp15_1 = A_01_2*w4;
2793                                      const register double tmp16_1 = tmp3_0*w0;                                      const double tmp16_1 = tmp3_0*w0;
2794                                      const register double tmp17_1 = A_01_1*w7;                                      const double tmp17_1 = A_01_1*w7;
2795                                      const register double tmp18_1 = tmp1_0*w6;                                      const double tmp18_1 = tmp1_0*w6;
2796                                      const register double tmp19_1 = A_10_3*w2;                                      const double tmp19_1 = A_10_3*w2;
2797                                      const register double tmp20_1 = tmp9_0*w4;                                      const double tmp20_1 = tmp9_0*w4;
2798                                      const register double tmp21_1 = tmp1_0*w11;                                      const double tmp21_1 = tmp1_0*w11;
2799                                      const register double tmp22_1 = tmp10_0*w7;                                      const double tmp22_1 = tmp10_0*w7;
2800                                      const register double tmp23_1 = tmp3_0*w14;                                      const double tmp23_1 = tmp3_0*w14;
2801                                      const register double tmp24_1 = tmp11_0*w13;                                      const double tmp24_1 = tmp11_0*w13;
2802                                      const register double tmp25_1 = tmp12_0*w12;                                      const double tmp25_1 = tmp12_0*w12;
2803                                      const register double tmp26_1 = tmp10_0*w4;                                      const double tmp26_1 = tmp10_0*w4;
2804                                      const register double tmp27_1 = tmp3_0*w11;                                      const double tmp27_1 = tmp3_0*w11;
2805                                      const register double tmp28_1 = tmp9_0*w7;                                      const double tmp28_1 = tmp9_0*w7;
2806                                      const register double tmp29_1 = tmp1_0*w14;                                      const double tmp29_1 = tmp1_0*w14;
2807                                      const register double tmp30_1 = tmp12_0*w13;                                      const double tmp30_1 = tmp12_0*w13;
2808                                      const register double tmp31_1 = tmp11_0*w12;                                      const double tmp31_1 = tmp11_0*w12;
2809                                      const register double tmp32_1 = tmp13_0*w2;                                      const double tmp32_1 = tmp13_0*w2;
2810                                      const register double tmp33_1 = tmp14_0*w5;                                      const double tmp33_1 = tmp14_0*w5;
2811                                      const register double tmp34_1 = tmp15_0*w8;                                      const double tmp34_1 = tmp15_0*w8;
2812                                      const register double tmp35_1 = A_01_0*w8;                                      const double tmp35_1 = A_01_0*w8;
2813                                      const register double tmp36_1 = tmp16_0*w1;                                      const double tmp36_1 = tmp16_0*w1;
2814                                      const register double tmp37_1 = A_10_1*w4;                                      const double tmp37_1 = A_10_1*w4;
2815                                      const register double tmp38_1 = tmp5_0*w15;                                      const double tmp38_1 = tmp5_0*w15;
2816                                      const register double tmp39_1 = A_10_2*w7;                                      const double tmp39_1 = A_10_2*w7;
2817                                      const register double tmp40_1 = tmp11_0*w17;                                      const double tmp40_1 = tmp11_0*w17;
2818                                      const register double tmp41_1 = A_01_3*w2;                                      const double tmp41_1 = A_01_3*w2;
2819                                      const register double tmp42_1 = tmp12_0*w16;                                      const double tmp42_1 = tmp12_0*w16;
2820                                      const register double tmp43_1 = tmp17_0*w5;                                      const double tmp43_1 = tmp17_0*w5;
2821                                      const register double tmp44_1 = tmp7_0*w7;                                      const double tmp44_1 = tmp7_0*w7;
2822                                      const register double tmp45_1 = tmp6_0*w4;                                      const double tmp45_1 = tmp6_0*w4;
2823                                      const register double tmp46_1 = A_01_3*w8;                                      const double tmp46_1 = A_01_3*w8;
2824                                      const register double tmp47_1 = A_10_2*w4;                                      const double tmp47_1 = A_10_2*w4;
2825                                      const register double tmp48_1 = A_10_1*w7;                                      const double tmp48_1 = A_10_1*w7;
2826                                      const register double tmp49_1 = tmp12_0*w17;                                      const double tmp49_1 = tmp12_0*w17;
2827                                      const register double tmp50_1 = A_01_0*w2;                                      const double tmp50_1 = A_01_0*w2;
2828                                      const register double tmp51_1 = tmp11_0*w16;                                      const double tmp51_1 = tmp11_0*w16;
2829                                      const register double tmp52_1 = tmp18_0*w8;                                      const double tmp52_1 = tmp18_0*w8;
2830                                      const register double tmp53_1 = tmp19_0*w2;                                      const double tmp53_1 = tmp19_0*w2;
2831                                      const register double tmp54_1 = tmp13_0*w8;                                      const double tmp54_1 = tmp13_0*w8;
2832                                      const register double tmp55_1 = tmp15_0*w2;                                      const double tmp55_1 = tmp15_0*w2;
2833                                      const register double tmp56_1 = tmp19_0*w8;                                      const double tmp56_1 = tmp19_0*w8;
2834                                      const register double tmp57_1 = tmp18_0*w2;                                      const double tmp57_1 = tmp18_0*w2;
2835                                      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,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
2836                                      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,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
2837                                      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,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
# Line 2607  void Rectangle::assemblePDESystem(Paso_S Line 2853  void Rectangle::assemblePDESystem(Paso_S
2853                          } else { // constant data                          } else { // constant data
2854                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2855                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2856                                      const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];                                      const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
2857                                      const register 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)];
2858                                      const register 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)];
2859                                      const register 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)];
2860                                      const register double tmp0_0 = A_01 + A_10;                                      const double tmp0_0 = A_01 + A_10;
2861                                      const register double tmp0_1 = A_00*w18;                                      const double tmp0_1 = A_00*w18;
2862                                      const register double tmp1_1 = A_01*w19;                                      const double tmp1_1 = A_01*w19;
2863                                      const register double tmp2_1 = A_10*w20;                                      const double tmp2_1 = A_10*w20;
2864                                      const register double tmp3_1 = A_11*w21;                                      const double tmp3_1 = A_11*w21;
2865                                      const register double tmp4_1 = A_00*w22;                                      const double tmp4_1 = A_00*w22;
2866                                      const register double tmp5_1 = tmp0_0*w19;                                      const double tmp5_1 = tmp0_0*w19;
2867                                      const register double tmp6_1 = A_11*w23;                                      const double tmp6_1 = A_11*w23;
2868                                      const register double tmp7_1 = A_11*w25;                                      const double tmp7_1 = A_11*w25;
2869                                      const register double tmp8_1 = A_00*w24;                                      const double tmp8_1 = A_00*w24;
2870                                      const register double tmp9_1 = tmp0_0*w20;                                      const double tmp9_1 = tmp0_0*w20;
2871                                      const register double tmp10_1 = A_01*w20;                                      const double tmp10_1 = A_01*w20;
2872                                      const register double tmp11_1 = A_11*w27;                                      const double tmp11_1 = A_11*w27;
2873                                      const register double tmp12_1 = A_00*w26;                                      const double tmp12_1 = A_00*w26;
2874                                      const register double tmp13_1 = A_10*w19;                                      const double tmp13_1 = A_10*w19;
2875                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2876                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2877                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
# Line 2655  void Rectangle::assemblePDESystem(Paso_S Line 2901  void Rectangle::assemblePDESystem(Paso_S
2901                          if (B.actsExpanded()) {                          if (B.actsExpanded()) {
2902                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2903                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2904                                      const register double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];                                      const double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
2905                                      const register double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];                                      const double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
2906                                      const register double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];                                      const double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
2907                                      const register double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];                                      const double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
2908                                      const register double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];                                      const double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];
2909                                      const register 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)];
2910                                      const register 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)];
2911                                      const register 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)];
2912                                      const register double tmp0_0 = B_1_0 + B_1_1;                                      const double tmp0_0 = B_1_0 + B_1_1;
2913                                      const register double tmp1_0 = B_1_2 + B_1_3;                                      const double tmp1_0 = B_1_2 + B_1_3;
2914                                      const register double tmp2_0 = B_0_1 + B_0_3;                                      const double tmp2_0 = B_0_1 + B_0_3;
2915                                      const register double tmp3_0 = B_0_0 + B_0_2;                                      const double tmp3_0 = B_0_0 + B_0_2;
2916                                      const register double tmp63_1 = B_1_1*w42;                                      const double tmp63_1 = B_1_1*w42;
2917                                      const register double tmp79_1 = B_1_1*w40;                                      const double tmp79_1 = B_1_1*w40;
2918                                      const register double tmp37_1 = tmp3_0*w35;                                      const double tmp37_1 = tmp3_0*w35;
2919                                      const register double tmp8_1 = tmp0_0*w32;                                      const double tmp8_1 = tmp0_0*w32;
2920                                      const register double tmp71_1 = B_0_1*w34;                                      const double tmp71_1 = B_0_1*w34;
2921                                      const register double tmp19_1 = B_0_3*w31;                                      const double tmp19_1 = B_0_3*w31;
2922                                      const register double tmp15_1 = B_0_3*w34;                                      const double tmp15_1 = B_0_3*w34;
2923                                      const register double tmp9_1 = tmp3_0*w34;                                      const double tmp9_1 = tmp3_0*w34;
2924                                      const register double tmp35_1 = B_1_0*w36;                                      const double tmp35_1 = B_1_0*w36;
2925                                      const register double tmp66_1 = B_0_3*w28;                                      const double tmp66_1 = B_0_3*w28;
2926                                      const register double tmp28_1 = B_1_0*w42;                                      const double tmp28_1 = B_1_0*w42;
2927                                      const register double tmp22_1 = B_1_0*w40;                                      const double tmp22_1 = B_1_0*w40;
2928                                      const register double tmp16_1 = B_1_2*w29;                                      const double tmp16_1 = B_1_2*w29;
2929                                      const register double tmp6_1 = tmp2_0*w35;                                      const double tmp6_1 = tmp2_0*w35;
2930                                      const register double tmp55_1 = B_1_3*w40;                                      const double tmp55_1 = B_1_3*w40;
2931                                      const register double tmp50_1 = B_1_3*w42;                                      const double tmp50_1 = B_1_3*w42;
2932                                      const register double tmp7_1 = tmp1_0*w29;                                      const double tmp7_1 = tmp1_0*w29;
2933                                      const register double tmp1_1 = tmp1_0*w32;                                      const double tmp1_1 = tmp1_0*w32;
2934                                      const register double tmp57_1 = B_0_3*w30;                                      const double tmp57_1 = B_0_3*w30;
2935                                      const register double tmp18_1 = B_1_1*w32;                                      const double tmp18_1 = B_1_1*w32;
2936                                      const register double tmp53_1 = B_1_0*w41;                                      const double tmp53_1 = B_1_0*w41;
2937                                      const register double tmp61_1 = B_1_3*w36;                                      const double tmp61_1 = B_1_3*w36;
2938                                      const register double tmp27_1 = B_0_3*w38;                                      const double tmp27_1 = B_0_3*w38;
2939                                      const register double tmp64_1 = B_0_2*w30;                                      const double tmp64_1 = B_0_2*w30;
2940                                      const register double tmp76_1 = B_0_1*w38;                                      const double tmp76_1 = B_0_1*w38;
2941                                      const register double tmp39_1 = tmp2_0*w34;                                      const double tmp39_1 = tmp2_0*w34;
2942                                      const register double tmp62_1 = B_0_1*w31;                                      const double tmp62_1 = B_0_1*w31;
2943                                      const register double tmp56_1 = B_0_0*w31;                                      const double tmp56_1 = B_0_0*w31;
2944                                      const register double tmp49_1 = B_1_1*w36;                                      const double tmp49_1 = B_1_1*w36;
2945                                      const register double tmp2_1 = B_0_2*w31;                                      const double tmp2_1 = B_0_2*w31;
2946                                      const register double tmp23_1 = B_0_2*w33;                                      const double tmp23_1 = B_0_2*w33;
2947                                      const register double tmp38_1 = B_1_1*w43;                                      const double tmp38_1 = B_1_1*w43;
2948                                      const register double tmp74_1 = B_1_2*w41;                                      const double tmp74_1 = B_1_2*w41;
2949                                      const register double tmp43_1 = B_1_1*w41;                                      const double tmp43_1 = B_1_1*w41;
2950                                      const register double tmp58_1 = B_0_2*w28;                                      const double tmp58_1 = B_0_2*w28;
2951                                      const register double tmp67_1 = B_0_0*w33;                                      const double tmp67_1 = B_0_0*w33;
2952                                      const register double tmp33_1 = tmp0_0*w39;                                      const double tmp33_1 = tmp0_0*w39;
2953                                      const register double tmp4_1 = B_0_0*w28;                                      const double tmp4_1 = B_0_0*w28;
2954                                      const register double tmp20_1 = B_0_0*w30;                                      const double tmp20_1 = B_0_0*w30;
2955                                      const register double tmp13_1 = B_0_2*w38;                                      const double tmp13_1 = B_0_2*w38;
2956                                      const register double tmp65_1 = B_1_2*w43;                                      const double tmp65_1 = B_1_2*w43;
2957                                      const register double tmp0_1 = tmp0_0*w29;                                      const double tmp0_1 = tmp0_0*w29;
2958                                      const register double tmp41_1 = tmp3_0*w33;                                      const double tmp41_1 = tmp3_0*w33;
2959                                      const register double tmp73_1 = B_0_2*w37;                                      const double tmp73_1 = B_0_2*w37;
2960                                      const register double tmp69_1 = B_0_0*w38;                                      const double tmp69_1 = B_0_0*w38;
2961                                      const register double tmp48_1 = B_1_2*w39;                                      const double tmp48_1 = B_1_2*w39;
2962                                      const register double tmp59_1 = B_0_1*w33;                                      const double tmp59_1 = B_0_1*w33;
2963                                      const register double tmp17_1 = B_1_3*w41;                                      const double tmp17_1 = B_1_3*w41;
2964                                      const register double tmp5_1 = B_0_3*w33;                                      const double tmp5_1 = B_0_3*w33;
2965                                      const register double tmp3_1 = B_0_1*w30;                                      const double tmp3_1 = B_0_1*w30;
2966                                      const register double tmp21_1 = B_0_1*w28;                                      const double tmp21_1 = B_0_1*w28;
2967                                      const register double tmp42_1 = B_1_0*w29;                                      const double tmp42_1 = B_1_0*w29;
2968                                      const register double tmp54_1 = B_1_2*w32;                                      const double tmp54_1 = B_1_2*w32;
2969                                      const register double tmp60_1 = B_1_0*w39;                                      const double tmp60_1 = B_1_0*w39;
2970                                      const register double tmp32_1 = tmp1_0*w36;                                      const double tmp32_1 = tmp1_0*w36;
2971                                      const register double tmp10_1 = B_0_1*w37;                                      const double tmp10_1 = B_0_1*w37;
2972                                      const register double tmp14_1 = B_0_0*w35;                                      const double tmp14_1 = B_0_0*w35;
2973                                      const register double tmp29_1 = B_0_1*w35;                                      const double tmp29_1 = B_0_1*w35;
2974                                      const register double tmp26_1 = B_1_2*w36;                                      const double tmp26_1 = B_1_2*w36;
2975                                      const register double tmp30_1 = B_1_3*w43;                                      const double tmp30_1 = B_1_3*w43;
2976                                      const register double tmp70_1 = B_0_2*w35;                                      const double tmp70_1 = B_0_2*w35;
2977                                      const register double tmp34_1 = B_1_3*w39;                                      const double tmp34_1 = B_1_3*w39;
2978                                      const register double tmp51_1 = B_1_0*w43;                                      const double tmp51_1 = B_1_0*w43;
2979                                      const register double tmp31_1 = B_0_2*w34;                                      const double tmp31_1 = B_0_2*w34;
2980                                      const register double tmp45_1 = tmp3_0*w28;                                      const double tmp45_1 = tmp3_0*w28;
2981                                      const register double tmp11_1 = tmp1_0*w39;                                      const double tmp11_1 = tmp1_0*w39;
2982                                      const register double tmp52_1 = B_1_1*w29;                                      const double tmp52_1 = B_1_1*w29;
2983                                      const register double tmp44_1 = B_1_3*w32;                                      const double tmp44_1 = B_1_3*w32;
2984                                      const register double tmp25_1 = B_1_1*w39;                                      const double tmp25_1 = B_1_1*w39;
2985                                      const register double tmp47_1 = tmp2_0*w33;                                      const double tmp47_1 = tmp2_0*w33;
2986                                      const register double tmp72_1 = B_1_3*w29;                                      const double tmp72_1 = B_1_3*w29;
2987                                      const register double tmp40_1 = tmp2_0*w28;                                      const double tmp40_1 = tmp2_0*w28;
2988                                      const register double tmp46_1 = B_1_2*w40;                                      const double tmp46_1 = B_1_2*w40;
2989                                      const register double tmp36_1 = B_1_2*w42;                                      const double tmp36_1 = B_1_2*w42;
2990                                      const register double tmp24_1 = B_0_0*w37;                                      const double tmp24_1 = B_0_0*w37;
2991                                      const register double tmp77_1 = B_0_3*w35;                                      const double tmp77_1 = B_0_3*w35;
2992                                      const register double tmp68_1 = B_0_3*w37;                                      const double tmp68_1 = B_0_3*w37;
2993                                      const register double tmp78_1 = B_0_0*w34;                                      const double tmp78_1 = B_0_0*w34;
2994                                      const register double tmp12_1 = tmp0_0*w36;                                      const double tmp12_1 = tmp0_0*w36;
2995                                      const register double tmp75_1 = B_1_0*w32;                                      const double tmp75_1 = B_1_0*w32;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
2996                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
2997                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2998                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2999                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
3000                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
3001                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
3002                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
3003                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
3004                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
3005                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
3006                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
3007                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
3008                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
3009                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
3010                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
3011                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
3012                                  }                                  }
3013                              }                              }
3014                          } else { // constant data                          } else { // constant data
3015                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3016                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3017                                      const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];                                      const double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
3018                                      const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];                                      const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
3019                                      const register double tmp6_1 = B_1*w50;                                      const double tmp6_1 = B_1*w50;
3020                                      const register double tmp1_1 = B_1*w45;                                      const double tmp1_1 = B_1*w45;
3021                                      const register double tmp5_1 = B_1*w49;                                      const double tmp5_1 = B_1*w49;
3022                                      const register double tmp4_1 = B_1*w48;                                      const double tmp4_1 = B_1*w48;
3023                                      const register double tmp0_1 = B_0*w44;                                      const double tmp0_1 = B_0*w44;
3024                                      const register double tmp2_1 = B_0*w46;                                      const double tmp2_1 = B_0*w46;
3025                                      const register double tmp7_1 = B_0*w51;                                      const double tmp7_1 = B_0*w51;
3026                                      const register double tmp3_1 = B_0*w47;                                      const double tmp3_1 = B_0*w47;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
3027                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;
3028                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3029                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;
3030                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
3031                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3032                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;
3033                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;
3034                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;
3035                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;
3036                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;  
3037                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;
3038                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
3039                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;
3040                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;
3041                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;
3042                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;
3043                                  }                                  }
3044                              }                              }
3045                          }                          }
# Line 2807  void Rectangle::assemblePDESystem(Paso_S Line 3053  void Rectangle::assemblePDESystem(Paso_S
3053                          if (C.actsExpanded()) {                          if (C.actsExpanded()) {
3054                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3055                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3056                                      const register double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];                                      const double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
3057                                      const register double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];                                      const double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
3058                                      const register double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];                                      const double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
3059                                      const register double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];                                      const double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
3060                                      const register double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];                                      const double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];
3061                                      const register double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];                                      const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
3062                                      const register double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];                                      const double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];
3063                                      const register double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];                                      const double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];
3064                                      const register double tmp2_0 = C_0_1 + C_0_3;                                      const double tmp2_0 = C_0_1 + C_0_3;
3065                                      const register double tmp1_0 = C_1_2 + C_1_3;                                      const double tmp1_0 = C_1_2 + C_1_3;
3066                                      const register double tmp3_0 = C_0_0 + C_0_2;                                      const double tmp3_0 = C_0_0 + C_0_2;
3067                                      const register double tmp0_0 = C_1_0 + C_1_1;                                      const double tmp0_0 = C_1_0 + C_1_1;
3068                                      const register double tmp64_1 = C_0_2*w30;                                      const double tmp64_1 = C_0_2*w30;
3069                                      const register double tmp14_1 = C_0_2*w28;                                      const double tmp14_1 = C_0_2*w28;
3070                                      const register double tmp19_1 = C_0_3*w31;                                      const double tmp19_1 = C_0_3*w31;
3071                                      const register double tmp22_1 = C_1_0*w40;                                      const double tmp22_1 = C_1_0*w40;
3072                                      const register double tmp37_1 = tmp3_0*w35;                                      const double tmp37_1 = tmp3_0*w35;
3073                                      const register double tmp29_1 = C_0_1*w35;                                      const double tmp29_1 = C_0_1*w35;
3074                                      const register double tmp73_1 = C_0_2*w37;                                      const double tmp73_1 = C_0_2*w37;
3075                                      const register double tmp74_1 = C_1_2*w41;                                      const double tmp74_1 = C_1_2*w41;
3076                                      const register double tmp52_1 = C_1_3*w39;                                      const double tmp52_1 = C_1_3*w39;
3077                                      const register double tmp25_1 = C_1_1*w39;                                      const double tmp25_1 = C_1_1*w39;
3078                                      const register double tmp62_1 = C_0_1*w31;                                      const double tmp62_1 = C_0_1*w31;
3079                                      const register double tmp79_1 = C_1_1*w40;                                      const double tmp79_1 = C_1_1*w40;
3080                                      const register double tmp43_1 = C_1_1*w36;                                      const double tmp43_1 = C_1_1*w36;
3081                                      const register double tmp27_1 = C_0_3*w38;                                      const double tmp27_1 = C_0_3*w38;
3082                                      const register double tmp28_1 = C_1_0*w42;                                      const double tmp28_1 = C_1_0*w42;
3083                                      const register double tmp63_1 = C_1_1*w42;                                      const double tmp63_1 = C_1_1*w42;
3084                                      const register double tmp59_1 = C_0_3*w34;                                      const double tmp59_1 = C_0_3*w34;
3085                                      const register double tmp72_1 = C_1_3*w29;                                      const double tmp72_1 = C_1_3*w29;
3086                                      const register double tmp40_1 = tmp2_0*w35;                                      const double tmp40_1 = tmp2_0*w35;
3087                                      const register double tmp13_1 = C_0_3*w30;                                      const double tmp13_1 = C_0_3*w30;
3088                                      const register double tmp51_1 = C_1_2*w40;                                      const double tmp51_1 = C_1_2*w40;
3089                                      const register double tmp54_1 = C_1_2*w42;                                      const double tmp54_1 = C_1_2*w42;
3090                                      const register double tmp12_1 = C_0_0*w31;                                      const double tmp12_1 = C_0_0*w31;
3091                                      const register double tmp2_1 = tmp1_0*w32;                                      const double tmp2_1 = tmp1_0*w32;
3092                                      const register double tmp68_1 = C_0_2*w31;                                      const double tmp68_1 = C_0_2*w31;
3093                                      const register double tmp75_1 = C_1_0*w32;                                      const double tmp75_1 = C_1_0*w32;
3094                                      const register double tmp49_1 = C_1_1*w41;                                      const double tmp49_1 = C_1_1*w41;
3095                                      const register double tmp4_1 = C_0_2*w35;                                      const double tmp4_1 = C_0_2*w35;
3096                                      const register double tmp66_1 = C_0_3*w28;                                      const double tmp66_1 = C_0_3*w28;
3097                                      const register double tmp56_1 = C_0_1*w37;                                      const double tmp56_1 = C_0_1*w37;
3098                                      const register double tmp5_1 = C_0_1*w34;                                      const double tmp5_1 = C_0_1*w34;
3099                                      const register double tmp38_1 = tmp2_0*w34;                                      const double tmp38_1 = tmp2_0*w34;
3100                                      const register double tmp76_1 = C_0_1*w38;                                      const double tmp76_1 = C_0_1*w38;
3101                                      const register double tmp21_1 = C_0_1*w28;                                      const double tmp21_1 = C_0_1*w28;
3102                                      const register double tmp69_1 = C_0_1*w30;                                      const double tmp69_1 = C_0_1*w30;
3103                                      const register double tmp53_1 = C_1_0*w36;                                      const double tmp53_1 = C_1_0*w36;
3104                                      const register double tmp42_1 = C_1_2*w39;                                      const double tmp42_1 = C_1_2*w39;
3105                                      const register double tmp32_1 = tmp1_0*w29;                                      const double tmp32_1 = tmp1_0*w29;
3106                                      const register double tmp45_1 = C_1_0*w43;                                      const double tmp45_1 = C_1_0*w43;
3107                                      const register double tmp33_1 = tmp0_0*w32;                                      const double tmp33_1 = tmp0_0*w32;
3108                                      const register double tmp35_1 = C_1_0*w41;                                      const double tmp35_1 = C_1_0*w41;
3109                                      const register double tmp26_1 = C_1_2*w36;                                      const double tmp26_1 = C_1_2*w36;
3110                                      const register double tmp67_1 = C_0_0*w33;                                      const double tmp67_1 = C_0_0*w33;
3111                                      const register double tmp31_1 = C_0_2*w34;                                      const double tmp31_1 = C_0_2*w34;
3112                                      const register double tmp20_1 = C_0_0*w30;                                      const double tmp20_1 = C_0_0*w30;
3113                                      const register double tmp70_1 = C_0_0*w28;                                      const double tmp70_1 = C_0_0*w28;
3114                                      const register double tmp8_1 = tmp0_0*w39;                                      const double tmp8_1 = tmp0_0*w39;
3115                                      const register double tmp30_1 = C_1_3*w43;                                      const double tmp30_1 = C_1_3*w43;
3116                                      const register double tmp0_1 = tmp0_0*w29;                                      const double tmp0_1 = tmp0_0*w29;
3117                                      const register double tmp17_1 = C_1_3*w41;                                      const double tmp17_1 = C_1_3*w41;
3118                                      const register double tmp58_1 = C_0_0*w35;                                      const double tmp58_1 = C_0_0*w35;
3119                                      const register double tmp9_1 = tmp3_0*w33;                                      const double tmp9_1 = tmp3_0*w33;
3120                                      const register double tmp61_1 = C_1_3*w36;                                      const double tmp61_1 = C_1_3*w36;
3121                                      const register double tmp41_1 = tmp3_0*w34;                                      const double tmp41_1 = tmp3_0*w34;
3122                                      const register double tmp50_1 = C_1_3*w32;                                      const double tmp50_1 = C_1_3*w32;
3123                                      const register double tmp18_1 = C_1_1*w32;                                      const double tmp18_1 = C_1_1*w32;
3124                                      const register double tmp6_1 = tmp1_0*w36;                                      const double tmp6_1 = tmp1_0*w36;
3125                                      const register double tmp3_1 = C_0_0*w38;                                      const double tmp3_1 = C_0_0*w38;
3126                                      const register double tmp34_1 = C_1_1*w29;                                      const double tmp34_1 = C_1_1*w29;
3127                                      const register double tmp77_1 = C_0_3*w35;                                      const double tmp77_1 = C_0_3*w35;
3128                                      const register double tmp65_1 = C_1_2*w43;                                      const double tmp65_1 = C_1_2*w43;
3129                                      const register double tmp71_1 = C_0_3*w33;                                      const double tmp71_1 = C_0_3*w33;
3130                                      const register double tmp55_1 = C_1_1*w43;                                      const double tmp55_1 = C_1_1*w43;
3131                                      const register double tmp46_1 = tmp3_0*w28;                                      const double tmp46_1 = tmp3_0*w28;
3132                                      const register double tmp24_1 = C_0_0*w37;                                      const double tmp24_1 = C_0_0*w37;
3133                                      const register double tmp10_1 = tmp1_0*w39;                                      const double tmp10_1 = tmp1_0*w39;
3134                                      const register double tmp48_1 = C_1_0*w29;                                      const double tmp48_1 = C_1_0*w29;
3135                                      const register double tmp15_1 = C_0_1*w33;                                      const double tmp15_1 = C_0_1*w33;
3136                                      const register double tmp36_1 = C_1_2*w32;                                      const double tmp36_1 = C_1_2*w32;
3137                                      const register double tmp60_1 = C_1_0*w39;                                      const double tmp60_1 = C_1_0*w39;
3138                                      const register double tmp47_1 = tmp2_0*w33;                                      const double tmp47_1 = tmp2_0*w33;
3139                                      const register double tmp16_1 = C_1_2*w29;                                      const double tmp16_1 = C_1_2*w29;
3140                                      const register double tmp1_1 = C_0_3*w37;                                      const double tmp1_1 = C_0_3*w37;
3141                                      const register double tmp7_1 = tmp2_0*w28;                                      const double tmp7_1 = tmp2_0*w28;
3142                                      const register double tmp39_1 = C_1_3*w40;                                      const double tmp39_1 = C_1_3*w40;
3143                                      const register double tmp44_1 = C_1_3*w42;                                      const double tmp44_1 = C_1_3*w42;
3144                                      const register double tmp57_1 = C_0_2*w38;                                      const double tmp57_1 = C_0_2*w38;
3145                                      const register double tmp78_1 = C_0_0*w34;                                      const double tmp78_1 = C_0_0*w34;
3146                                      const register double tmp11_1 = tmp0_0*w36;                                      const double tmp11_1 = tmp0_0*w36;
3147                                      const register double tmp23_1 = C_0_2*w33;                                      const double tmp23_1 = C_0_2*w33;
3148                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
3149                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
3150                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
# Line 2920  void Rectangle::assemblePDESystem(Paso_S Line 3166  void Rectangle::assemblePDESystem(Paso_S
3166                          } else { // constant data                          } else { // constant data
3167                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3168                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3169                                      const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
3170                                      const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
3171                                      const register double tmp1_1 = C_1*w45;                                      const double tmp1_1 = C_1*w45;
3172                                      const register double tmp3_1 = C_0*w51;                                      const double tmp3_1 = C_0*w51;
3173                                      const register double tmp4_1 = C_0*w44;                                      const double tmp4_1 = C_0*w44;
3174                                      const register double tmp7_1 = C_0*w46;                                      const double tmp7_1 = C_0*w46;
3175                                      const register double tmp5_1 = C_1*w49;                                      const double tmp5_1 = C_1*w49;
3176                                      const register double tmp2_1 = C_1*w48;                                      const double tmp2_1 = C_1*w48;
3177                                      const register double tmp0_1 = C_0*w47;                                      const double tmp0_1 = C_0*w47;
3178                                      const register double tmp6_1 = C_1*w50;                                      const double tmp6_1 = C_1*w50;
3179                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3180                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3181                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
# Line 2959  void Rectangle::assemblePDESystem(Paso_S Line 3205  void Rectangle::assemblePDESystem(Paso_S
3205                          if (D.actsExpanded()) {                          if (D.actsExpanded()) {
3206                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3207                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3208                                      const register double D_0 = D_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double D_0 = D_p[INDEX3(k, m, 0, numEq, numComp)];
3209                                      const register double D_1 = D_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double D_1 = D_p[INDEX3(k, m, 1, numEq, numComp)];
3210                                      const register double D_2 = D_p[INDEX3(k, m, 2, numEq, numComp)];                                      const double D_2 = D_p[INDEX3(k, m, 2, numEq, numComp)];
3211                                      const register double D_3 = D_p[INDEX3(k, m, 3, numEq, numComp)];                                      const double D_3 = D_p[INDEX3(k, m, 3, numEq, numComp)];
3212                                      const register double tmp4_0 = D_1 + D_3;                                      const double tmp4_0 = D_1 + D_3;
3213                                      const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;                                      const double tmp2_0 = D_0 + D_1 + D_2 + D_3;
3214                                      const register double tmp5_0 = D_0 + D_2;                                      const double tmp5_0 = D_0 + D_2;
3215                                      const register double tmp0_0 = D_0 + D_1;                                      const double tmp0_0 = D_0 + D_1;
3216                                      const register double tmp6_0 = D_0 + D_3;                                      const double tmp6_0 = D_0 + D_3;
3217                                      const register double tmp1_0 = D_2 + D_3;                                      const double tmp1_0 = D_2 + D_3;
3218                                      const register double tmp3_0 = D_1 + D_2;                                      const double tmp3_0 = D_1 + D_2;
3219                                      const register double tmp16_1 = D_1*w56;                                      const double tmp16_1 = D_1*w56;
3220                                      const register double tmp14_1 = tmp6_0*w54;                                      const double tmp14_1 = tmp6_0*w54;
3221                                      const register double tmp8_1 = D_3*w55;                                      const double tmp8_1 = D_3*w55;
3222                                      const register double tmp2_1 = tmp2_0*w54;                                      const double tmp2_1 = tmp2_0*w54;
3223                                      const register double tmp12_1 = tmp5_0*w52;                                      const double tmp12_1 = tmp5_0*w52;
3224                                      const register double tmp4_1 = tmp0_0*w53;                                      const double tmp4_1 = tmp0_0*w53;
3225                                      const register double tmp3_1 = tmp1_0*w52;                                      const double tmp3_1 = tmp1_0*w52;
3226                                      const register double tmp13_1 = tmp4_0*w53;                                      const double tmp13_1 = tmp4_0*w53;
3227                                      const register double tmp10_1 = tmp4_0*w52;                                      const double tmp10_1 = tmp4_0*w52;
3228                                      const register double tmp1_1 = tmp1_0*w53;                                      const double tmp1_1 = tmp1_0*w53;
3229                                      const register double tmp7_1 = D_3*w56;                                      const double tmp7_1 = D_3*w56;
3230                                      const register double tmp0_1 = tmp0_0*w52;                                      const double tmp0_1 = tmp0_0*w52;
3231                                      const register double tmp11_1 = tmp5_0*w53;                                      const double tmp11_1 = tmp5_0*w53;
3232                                      const register double tmp9_1 = D_0*w56;                                      const double tmp9_1 = D_0*w56;
3233                                      const register double tmp5_1 = tmp3_0*w54;                                      const double tmp5_1 = tmp3_0*w54;
3234                                      const register double tmp18_1 = D_2*w56;                                      const double tmp18_1 = D_2*w56;
3235                                      const register double tmp17_1 = D_1*w55;                                      const double tmp17_1 = D_1*w55;
3236                                      const register double tmp6_1 = D_0*w55;                                      const double tmp6_1 = D_0*w55;
3237                                      const register double tmp15_1 = D_2*w55;                                      const double tmp15_1 = D_2*w55;
3238                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3239                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1;
3240                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
# Line 3010  void Rectangle::assemblePDESystem(Paso_S Line 3256  void Rectangle::assemblePDESystem(Paso_S
3256                          } else { // constant data                          } else { // constant data
3257                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3258                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3259                                      const register double D_0 = D_p[INDEX2(k, m, numEq)];                                      const double D_0 = D_p[INDEX2(k, m, numEq)];
3260                                      const register double tmp2_1 = D_0*w59;                                      const double tmp2_1 = D_0*w59;
3261                                      const register double tmp1_1 = D_0*w58;                                      const double tmp1_1 = D_0*w58;
3262                                      const register double tmp0_1 = D_0*w57;                                      const double tmp0_1 = D_0*w57;
3263                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
3264                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1;
3265                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
# Line 3042  void Rectangle::assemblePDESystem(Paso_S Line 3288  void Rectangle::assemblePDESystem(Paso_S
3288                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3289                          if (X.actsExpanded()) {                          if (X.actsExpanded()) {
3290                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3291                                  const register double X_0_0 = X_p[INDEX3(k, 0, 0, numEq, 2)];                                  const double X_0_0 = X_p[INDEX3(k, 0, 0, numEq, 2)];
3292                                  const register double X_1_0 = X_p[INDEX3(k, 1, 0, numEq, 2)];                                  const double X_1_0 = X_p[INDEX3(k, 1, 0, numEq, 2)];
3293                                  const register double X_0_1 = X_p[INDEX3(k, 0, 1, numEq, 2)];                                  const double X_0_1 = X_p[INDEX3(k, 0, 1, numEq, 2)];
3294                                  const register double X_1_1 = X_p[INDEX3(k, 1, 1, numEq, 2)];                                  const double X_1_1 = X_p[INDEX3(k, 1, 1, numEq, 2)];
3295                                  const register double X_0_2 = X_p[INDEX3(k, 0, 2, numEq, 2)];                                  const double X_0_2 = X_p[INDEX3(k, 0, 2, numEq, 2)];
3296                                  const register double X_1_2 = X_p[INDEX3(k, 1, 2, numEq, 2)];                                  const double X_1_2 = X_p[INDEX3(k, 1, 2, numEq, 2)];
3297                                  const register double X_0_3 = X_p[INDEX3(k, 0, 3, numEq, 2)];                                  const double X_0_3 = X_p[INDEX3(k, 0, 3, numEq, 2)];
3298                                  const register double X_1_3 = X_p[INDEX3(k, 1, 3, numEq, 2)];                                  const double X_1_3 = X_p[INDEX3(k, 1, 3, numEq, 2)];
3299                                  const register double tmp1_0 = X_1_1 + X_1_3;                                  const double tmp1_0 = X_1_1 + X_1_3;
3300                                  const register double tmp3_0 = X_0_0 + X_0_1;                                  const double tmp3_0 = X_0_0 + X_0_1;
3301                                  const register double tmp2_0 = X_1_0 + X_1_2;                                  const double tmp2_0 = X_1_0 + X_1_2;
3302                                  const register double tmp0_0 = X_0_2 + X_0_3;                                  const double tmp0_0 = X_0_2 + X_0_3;
3303                                  const register double tmp8_1 = tmp2_0*w66;                                  const double tmp8_1 = tmp2_0*w66;
3304                                  const register double tmp5_1 = tmp3_0*w64;                                  const double tmp5_1 = tmp3_0*w64;
3305                                  const register double tmp14_1 = tmp0_0*w64;                                  const double tmp14_1 = tmp0_0*w64;
3306                                  const register double tmp3_1 = tmp3_0*w60;                                  const double tmp3_1 = tmp3_0*w60;
3307                                  const register double tmp9_1 = tmp3_0*w63;                                  const double tmp9_1 = tmp3_0*w63;
3308                                  const register double tmp13_1 = tmp3_0*w65;                                  const double tmp13_1 = tmp3_0*w65;
3309                                  const register double tmp12_1 = tmp1_0*w66;                                  const double tmp12_1 = tmp1_0*w66;
3310                                  const register double tmp10_1 = tmp0_0*w60;                                  const double tmp10_1 = tmp0_0*w60;
3311                                  const register double tmp2_1 = tmp2_0*w61;                                  const double tmp2_1 = tmp2_0*w61;
3312                                  const register double tmp6_1 = tmp2_0*w62;                                  const double tmp6_1 = tmp2_0*w62;
3313                                  const register double tmp4_1 = tmp0_0*w65;                                  const double tmp4_1 = tmp0_0*w65;
3314                                  const register double tmp11_1 = tmp1_0*w67;                                  const double tmp11_1 = tmp1_0*w67;
3315                                  const register double tmp1_1 = tmp1_0*w62;                                  const double tmp1_1 = tmp1_0*w62;
3316                                  const register double tmp7_1 = tmp1_0*w61;                                  const double tmp7_1 = tmp1_0*w61;
3317                                  const register double tmp0_1 = tmp0_0*w63;                                  const double tmp0_1 = tmp0_0*w63;
3318                                  const register double tmp15_1 = tmp2_0*w67;                                  const double tmp15_1 = tmp2_0*w67;
3319                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3320                                  EM_F[INDEX2(k,1,numEq)]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
3321                                  EM_F[INDEX2(k,2,numEq)]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
# Line 3077  void Rectangle::assemblePDESystem(Paso_S Line 3323  void Rectangle::assemblePDESystem(Paso_S
3323                              }                              }
3324                          } else { // constant data                          } else { // constant data
3325                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3326                                  const register double X_0 = X_p[INDEX2(k, 0, numEq)];                                  const double X_0 = X_p[INDEX2(k, 0, numEq)];
3327                                  const register double X_1 = X_p[INDEX2(k, 1, numEq)];                                  const double X_1 = X_p[INDEX2(k, 1, numEq)];
3328                                  const register double tmp0_1 = X_1*w69;                                  const double tmp0_1 = X_1*w69;
3329                                  const register double tmp1_1 = X_0*w68;                                  const double tmp1_1 = X_0*w68;
3330                                  const register double tmp2_1 = X_0*w70;                                  const double tmp2_1 = X_0*w70;
3331                                  const register double tmp3_1 = X_1*w71;                                  const double tmp3_1 = X_1*w71;
3332                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3333                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;
3334                                  EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;
# Line 3098  void Rectangle::assemblePDESystem(Paso_S Line 3344  void Rectangle::assemblePDESystem(Paso_S
3344                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3345                          if (Y.actsExpanded()) {                          if (Y.actsExpanded()) {
3346                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3347                                  const register double Y_0 = Y_p[INDEX2(k, 0, numEq)];                                  const double Y_0 = Y_p[INDEX2(k, 0, numEq)];
3348                                  const register double Y_1 = Y_p[INDEX2(k, 1, numEq)];                                  const double Y_1 = Y_p[INDEX2(k, 1, numEq)];
3349                                  const register double Y_2 = Y_p[INDEX2(k, 2, numEq)];                                  const double Y_2 = Y_p[INDEX2(k, 2, numEq)];
3350                                  const register double Y_3 = Y_p[INDEX2(k, 3, numEq)];                                  const double Y_3 = Y_p[INDEX2(k, 3, numEq)];
3351                                  const register double tmp0_0 = Y_1 + Y_2;                                  const double tmp0_0 = Y_1 + Y_2;
3352                                  const register double tmp1_0 = Y_0 + Y_3;                                  const double tmp1_0 = Y_0 + Y_3;
3353                                  const register double tmp0_1 = Y_0*w72;                                  const double tmp0_1 = Y_0*w72;
3354                                  const register double tmp1_1 = tmp0_0*w73;                                  const double tmp1_1 = tmp0_0*w73;
3355                                  const register double tmp2_1 = Y_3*w74;                                  const double tmp2_1 = Y_3*w74;
3356                                  const register double tmp3_1 = Y_1*w72;                                  const double tmp3_1 = Y_1*w72;
3357                                  const register double tmp4_1 = tmp1_0*w73;                                  const double tmp4_1 = tmp1_0*w73;
3358                                  const register double tmp5_1 = Y_2*w74;                                  const double tmp5_1 = Y_2*w74;
3359                                  const register double tmp6_1 = Y_2*w72;                                  const double tmp6_1 = Y_2*w72;
3360                                  const register double tmp7_1 = Y_1*w74;                                  const double tmp7_1 = Y_1*w74;
3361                                  const register double tmp8_1 = Y_3*w72;                                  const double tmp8_1 = Y_3*w72;
3362                                  const register double tmp9_1 = Y_0*w74;                                  const double tmp9_1 = Y_0*w74;
3363                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
3364                                  EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
3365                                  EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
# Line 3121  void Rectangle::assemblePDESystem(Paso_S Line 3367  void Rectangle::assemblePDESystem(Paso_S
3367                              }                              }
3368                          } else { // constant data                          } else { // constant data
3369                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3370                                  const register double tmp0_1 = Y_p[k]*w75;                                  const double tmp0_1 = Y_p[k]*w75;
3371                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3372                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3373                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
# Line 3132  void Rectangle::assemblePDESystem(Paso_S Line 3378  void Rectangle::assemblePDESystem(Paso_S
3378    
3379                      // 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)
3380                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
3381                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
3382                      rowIndex.push_back(m_dofMap[firstNode]);                              firstNode, numEq, numComp);
                     rowIndex.push_back(m_dofMap[firstNode+1]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                     if (add_EM_F) {  
                         double *F_p=rhs.getSampleDataRW(0);  
                         for (index_t i=0; i<rowIndex.size(); i++) {  
                             if (rowIndex[i]<getNumDOF()) {  
                                 for (index_t eq=0; eq<numEq; eq++) {  
                                     F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];  
                                 }  
                             }  
                         }  
                     }  
                     if (add_EM_S) {  
                         addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);  
                     }  
   
3383                  } // end k0 loop                  } // end k0 loop
3384              } // end k1 loop              } // end k1 loop
3385          } // end of colouring          } // end of colouring
# Line 3161  void Rectangle::assemblePDESystem(Paso_S Line 3390  void Rectangle::assemblePDESystem(Paso_S
3390  void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,
3391          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
3392          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
3393          const escript::Data& X, const escript::Data& Y,          const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
3394  {  {
3395      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
3396      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 3211  void Rectangle::assemblePDESystemReduced Line 3439  void Rectangle::assemblePDESystemReduced
3439                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
3440                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3441                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3442                                  const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];                                  const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
3443                                  const register 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)];
3444                                  const register 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)];
3445                                  const register 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)];
3446                                  const register double tmp0_0 = A_01 + A_10;                                  const double tmp0_0 = A_01 + A_10;
3447                                  const register double tmp0_1 = A_11*w3;                                  const double tmp0_1 = A_11*w3;
3448                                  const register double tmp1_1 = A_00*w0;                                  const double tmp1_1 = A_00*w0;
3449                                  const register double tmp2_1 = A_01*w1;                                  const double tmp2_1 = A_01*w1;
3450                                  const register double tmp3_1 = A_10*w2;                                  const double tmp3_1 = A_10*w2;
3451                                  const register double tmp4_1 = tmp0_0*w1;                                  const double tmp4_1 = tmp0_0*w1;
3452                                  const register double tmp5_1 = A_11*w4;                                  const double tmp5_1 = A_11*w4;
3453                                  const register double tmp6_1 = A_00*w5;                                  const double tmp6_1 = A_00*w5;
3454                                  const register double tmp7_1 = tmp0_0*w2;                                  const double tmp7_1 = tmp0_0*w2;
3455                                  const register double tmp8_1 = A_10*w1;                                  const double tmp8_1 = A_10*w1;
3456                                  const register double tmp9_1 = A_01*w2;                                  const double tmp9_1 = A_01*w2;
3457                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
3458                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
3459                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
# Line 3253  void Rectangle::assemblePDESystemReduced Line 3481  void Rectangle::assemblePDESystemReduced
3481                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
3482                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3483                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3484                                  const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];                                  const double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
3485                                  const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];                                  const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
3486                                  const register double tmp0_1 = B_0*w6;                                  const double tmp0_1 = B_0*w6;
3487                                  const register double tmp1_1 = B_1*w7;                                  const double tmp1_1 = B_1*w7;
3488                                  const register double tmp2_1 = B_0*w8;                                  const double tmp2_1 = B_0*w8;
3489                                  const register double tmp3_1 = B_1*w9;                                  const double tmp3_1 = B_1*w9;
3490                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3491                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3492                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
# Line 3286  void Rectangle::assemblePDESystemReduced Line 3514  void Rectangle::assemblePDESystemReduced
3514                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
3515                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3516                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3517                                  const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];                                  const double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
3518                                  const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];                                  const double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
3519                                  const register double tmp0_1 = C_0*w8;                                  const double tmp0_1 = C_0*w8;
3520                                  const register double tmp1_1 = C_1*w7;                                  const double tmp1_1 = C_1*w7;
3521                                  const register double tmp2_1 = C_1*w9;                                  const double tmp2_1 = C_1*w9;
3522                                  const register double tmp3_1 = C_0*w6;                                  const double tmp3_1 = C_0*w6;
3523                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3524                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3525                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
# Line 3319  void Rectangle::assemblePDESystemReduced Line 3547  void Rectangle::assemblePDESystemReduced
3547                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
3548                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3549                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3550                                  const register double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;                                  const double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;
3551                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
3552                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
3553                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;
# Line 3346  void Rectangle::assemblePDESystemReduced Line 3574  void Rectangle::assemblePDESystemReduced
3574                          add_EM_F=true;                          add_EM_F=true;
3575                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3576                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3577                              const register double X_0 = X_p[INDEX2(k, 0, numEq)];                              const double X_0 = X_p[INDEX2(k, 0, numEq)];
3578                              const register double X_1 = X_p[INDEX2(k, 1, numEq)];                              const double X_1 = X_p[INDEX2(k, 1, numEq)];
3579                              const register double tmp0_1 = X_0*w11;                              const double tmp0_1 = X_0*w11;
3580                              const register double tmp1_1 = X_1*w12;                              const double tm