/[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 3744 by caltinay, Tue Dec 13 06:41:54 2011 UTC trunk/ripley/src/Rectangle.cpp revision 4002 by caltinay, Fri Sep 28 00:16:56 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2011 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/SystemMatrixPattern.h"  #include <paso/SystemMatrix.h>
19  }  }
20    
21  #if USE_SILO  #if USE_SILO
# Line 29  using namespace std; Line 31  using namespace std;
31    
32  namespace ripley {  namespace ripley {
33    
34  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,
35                         double y1, int d0, int d1) :
36      RipleyDomain(2),      RipleyDomain(2),
37      m_gNE0(n0),      m_x0(x0),
38      m_gNE1(n1),      m_y0(y0),
39      m_l0(l0),      m_l0(x1-x0),
40      m_l1(l1),      m_l1(y1-y0)
     m_NX(d0),  
     m_NY(d1)  
41  {  {
42        // ignore subdivision parameters for serial run
43        if (m_mpiInfo->size == 1) {
44            d0=1;
45            d1=1;
46        }
47    
48        bool warn=false;
49        // if number of subdivisions is non-positive, try to subdivide by the same
50        // ratio as the number of elements
51        if (d0<=0 && d1<=0) {
52            warn=true;
53            d0=(int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1)));
54            d1=m_mpiInfo->size/d0;
55            if (d0*d1 != m_mpiInfo->size) {
56                // ratios not the same so subdivide side with more elements only
57                if (n0>n1) {
58                    d0=0;
59                    d1=1;
60                } else {
61                    d0=1;
62                    d1=0;
63                }
64            }
65        }
66        if (d0<=0) {
67            // d1 is preset, determine d0 - throw further down if result is no good
68            d0=m_mpiInfo->size/d1;
69        } else if (d1<=0) {
70            // d0 is preset, determine d1 - throw further down if result is no good
71            d1=m_mpiInfo->size/d0;
72        }
73    
74        m_NX=d0;
75        m_NY=d1;
76    
77      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
78      // among number of ranks      // among number of ranks
79      if (m_NX*m_NY != m_mpiInfo->size)      if (m_NX*m_NY != m_mpiInfo->size)
80          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
81    
82      if (n0%m_NX > 0 || n1%m_NY > 0)      if (warn) {
83          throw RipleyException("Number of elements must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
84                << d1 << "). This may not be optimal!" << endl;
85        }
86    
87        if ((n0+1)%m_NX > 0) {
88            double Dx=m_l0/n0;
89            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
90            m_l0=Dx*n0;
91            cout << "Warning: Adjusted number of elements and length. N0="
92                << n0 << ", l0=" << m_l0 << endl;
93        }
94        if ((n1+1)%m_NY > 0) {
95            double Dy=m_l1/n1;
96            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
97            m_l1=Dy*n1;
98            cout << "Warning: Adjusted number of elements and length. N1="
99                << n1 << ", l1=" << m_l1 << endl;
100        }
101    
102        m_gNE0=n0;
103        m_gNE1=n1;
104    
105        if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
106            throw RipleyException("Too few elements for the number of ranks");
107    
108        // local number of elements (with and without overlap)
109        m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
110        if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
111            m_NE0++;
112        else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)
113            m_ownNE0--;
114    
115        m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
116        if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)
117            m_NE1++;
118        else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)
119            m_ownNE1--;
120    
121      // local number of elements      // local number of nodes
     m_NE0 = n0/m_NX;  
     m_NE1 = n1/m_NY;  
     // local number of nodes (not necessarily owned)  
122      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
123      m_N1 = m_NE1+1;      m_N1 = m_NE1+1;
124    
125      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
126      m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);
127      m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);      if (m_offset0 > 0)
128            m_offset0--;
129        m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);
130        if (m_offset1 > 0)
131            m_offset1--;
132    
133      populateSampleIds();      populateSampleIds();
134        createPattern();
135  }  }
136    
137  Rectangle::~Rectangle()  Rectangle::~Rectangle()
138  {  {
139        Paso_SystemMatrixPattern_free(m_pattern);
140        Paso_Connector_free(m_connector);
141  }  }
142    
143  string Rectangle::getDescription() const  string Rectangle::getDescription() const
# Line 73  bool Rectangle::operator==(const Abstrac Line 151  bool Rectangle::operator==(const Abstrac
151      if (o) {      if (o) {
152          return (RipleyDomain::operator==(other) &&          return (RipleyDomain::operator==(other) &&
153                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1
154                    && m_x0==o->m_x0 && m_y0==o->m_y0
155                  && m_l0==o->m_l0 && m_l1==o->m_l1                  && m_l0==o->m_l0 && m_l1==o->m_l1
156                  && m_NX==o->m_NX && m_NY==o->m_NY);                  && m_NX==o->m_NX && m_NY==o->m_NY);
157      }      }
# Line 80  bool Rectangle::operator==(const Abstrac Line 159  bool Rectangle::operator==(const Abstrac
159      return false;      return false;
160  }  }
161    
162    void Rectangle::readBinaryGrid(escript::Data& out, string filename,
163                const vector<int>& first, const vector<int>& numValues) const
164    {
165        // check destination function space
166        int myN0, myN1;
167        if (out.getFunctionSpace().getTypeCode() == Nodes) {
168            myN0 = m_N0;
169            myN1 = m_N1;
170        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
171                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
172            myN0 = m_NE0;
173            myN1 = m_NE1;
174        } else
175            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
176    
177        // check file existence and size
178        ifstream f(filename.c_str(), ifstream::binary);
179        if (f.fail()) {
180            throw RipleyException("readBinaryGrid(): cannot open file");
181        }
182        f.seekg(0, ios::end);
183        const int numComp = out.getDataPointSize();
184        const int filesize = f.tellg();
185        const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(float);
186        if (filesize < reqsize) {
187            f.close();
188            throw RipleyException("readBinaryGrid(): not enough data in file");
189        }
190    
191        // check if this rank contributes anything
192        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
193                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1) {
194            f.close();
195            return;
196        }
197    
198        // now determine how much this rank has to write
199    
200        // first coordinates in data object to write to
201        const int first0 = max(0, first[0]-m_offset0);
202        const int first1 = max(0, first[1]-m_offset1);
203        // indices to first value in file
204        const int idx0 = max(0, m_offset0-first[0]);
205        const int idx1 = max(0, m_offset1-first[1]);
206        // number of values to write
207        const int num0 = min(numValues[0]-idx0, myN0-first0);
208        const int num1 = min(numValues[1]-idx1, myN1-first1);
209    
210        out.requireWrite();
211        vector<float> values(num0*numComp);
212        const int dpp = out.getNumDataPointsPerSample();
213    
214        for (index_t y=0; y<num1; y++) {
215            const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);
216            f.seekg(fileofs*sizeof(float));
217            f.read((char*)&values[0], num0*numComp*sizeof(float));
218            for (index_t x=0; x<num0; x++) {
219                double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0);
220                for (index_t c=0; c<numComp; c++) {
221                    for (index_t q=0; q<dpp; q++) {
222                        *dest++ = static_cast<double>(values[x*numComp+c]);
223                    }
224                }
225            }
226        }
227    
228        f.close();
229    }
230    
231  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
232  {  {
233  #if USE_SILO  #if USE_SILO
# Line 88  void Rectangle::dump(const string& fileN Line 236  void Rectangle::dump(const string& fileN
236          fn+=".silo";          fn+=".silo";
237      }      }
238    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
239      int driver=DB_HDF5;          int driver=DB_HDF5;    
     string siloPath;  
240      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
241        const char* blockDirFmt = "/block%04d";
242    
243  #ifdef ESYS_MPI  #ifdef ESYS_MPI
244      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
245        const int NUM_SILO_FILES = 1;
246  #endif  #endif
247    
248      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 111  void Rectangle::dump(const string& fileN Line 258  void Rectangle::dump(const string& fileN
258                          PMPIO_DefaultClose, (void*)&driver);                          PMPIO_DefaultClose, (void*)&driver);
259          }          }
260          if (baton) {          if (baton) {
261              char str[64];              char siloPath[64];
262              snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));              snprintf(siloPath, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
263              siloPath = str;              dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath);
             dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());  
264          }          }
265  #endif  #endif
266      } else {      } else {
# Line 126  void Rectangle::dump(const string& fileN Line 272  void Rectangle::dump(const string& fileN
272              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
273                      getDescription().c_str(), driver);                      getDescription().c_str(), driver);
274          }          }
275            char siloPath[64];
276            snprintf(siloPath, 64, blockDirFmt, 0);
277            DBMkDir(dbfile, siloPath);
278            DBSetDir(dbfile, siloPath);
279      }      }
280    
281      if (!dbfile)      if (!dbfile)
# Line 218  void Rectangle::dump(const string& fileN Line 368  void Rectangle::dump(const string& fileN
368      }      }
369    
370  #else // USE_SILO  #else // USE_SILO
371      throw RipleyException("dump(): no Silo support");      throw RipleyException("dump: no Silo support");
372  #endif  #endif
373  }  }
374    
# Line 226  const int* Rectangle::borrowSampleRefere Line 376  const int* Rectangle::borrowSampleRefere
376  {  {
377      switch (fsType) {      switch (fsType) {
378          case Nodes:          case Nodes:
379            case ReducedNodes: // FIXME: reduced
380              return &m_nodeId[0];              return &m_nodeId[0];
381          case DegreesOfFreedom:          case DegreesOfFreedom:
382          case ReducedDegreesOfFreedom: //FIXME: reduced          case ReducedDegreesOfFreedom: // FIXME: reduced
383              return &m_dofId[0];              return &m_dofId[0];
384          case Elements:          case Elements:
385          case ReducedElements:          case ReducedElements:
# Line 241  const int* Rectangle::borrowSampleRefere Line 392  const int* Rectangle::borrowSampleRefere
392      }      }
393    
394      stringstream msg;      stringstream msg;
395      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
396      throw RipleyException(msg.str());      throw RipleyException(msg.str());
397  }  }
398    
399  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsType, index_t id) const
400  {  {
401  #ifdef ESYS_MPI      if (getMPISize()==1)
402      if (fsCode != ReducedNodes) {          return true;
         const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];  
         const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;  
         return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);  
     } else {  
         stringstream msg;  
         msg << "ownSample() not implemented for "  
             << functionSpaceTypeAsString(fsCode);  
         throw RipleyException(msg.str());  
     }  
 #else  
     return true;  
 #endif  
 }  
   
 void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const  
 {  
     escript::Data& in = *const_cast<escript::Data*>(&cIn);  
     const dim_t numComp = in.getDataPointSize();  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     const double cx0 = -1./h0;  
     const double cx1 = -.78867513459481288225/h0;  
     const double cx2 = -.5/h0;  
     const double cx3 = -.21132486540518711775/h0;  
     const double cx4 = .21132486540518711775/h0;  
     const double cx5 = .5/h0;  
     const double cx6 = .78867513459481288225/h0;  
     const double cx7 = 1./h0;  
     const double cy0 = -1./h1;  
     const double cy1 = -.78867513459481288225/h1;  
     const double cy2 = -.5/h1;  
     const double cy3 = -.21132486540518711775/h1;  
     const double cy4 = .21132486540518711775/h1;  
     const double cy5 = .5/h1;  
     const double cy6 = .78867513459481288225/h1;  
     const double cy7 = 1./h1;  
403    
404      if (out.getFunctionSpace().getTypeCode() == Elements) {      switch (fsType) {
405          /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */          case Nodes:
406  #pragma omp parallel for          case ReducedNodes: // FIXME: reduced
407          for (index_t k1=0; k1 < m_NE1; ++k1) {              return (m_dofMap[id] < getNumDOF());
408              for (index_t k0=0; k0 < m_NE0; ++k0) {          case DegreesOfFreedom:
409                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));          case ReducedDegreesOfFreedom:
410                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              return true;
411                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));          case Elements:
412                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));          case ReducedElements:
413                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              // check ownership of element's bottom left node
414                  for (index_t i=0; i < numComp; ++i) {              return (m_dofMap[id%m_NE0+m_N0*(id/m_NE0)] < getNumDOF());
415                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;          case FaceElements:
416                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;          case ReducedFaceElements:
417                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;              {
418                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                  // determine which face the sample belongs to before
419                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                  // checking ownership of corresponding element's first node
420                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                  const IndexVector faces = getNumFacesPerBoundary();
421                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                  dim_t n=0;
422                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                  for (size_t i=0; i<faces.size(); i++) {
423                  } /* end of component loop i */                      n+=faces[i];
424              } /* end of k0 loop */                      if (id<n) {
425          } /* end of k1 loop */                          index_t k;
426          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */                          if (i==1)
427      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {                              k=m_N0-2;
428          /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */                          else if (i==3)
429  #pragma omp parallel for                              k=m_N0*(m_N1-2);
430          for (index_t k1=0; k1 < m_NE1; ++k1) {                          else
431              for (index_t k0=0; k0 < m_NE0; ++k0) {                              k=0;
432                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                          // determine whether to move right or up
433                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                          const index_t delta=(i/2==0 ? m_N0 : 1);
434                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                          return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());
435                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                      }
436                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  }
437                  for (index_t i=0; i < numComp; ++i) {                  return false;
                     o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);  
                     o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);  
                 } /* end of component loop i */  
             } /* end of k0 loop */  
         } /* end of k1 loop */  
         /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */  
     } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {  
 #pragma omp parallel  
         {  
             /*** GENERATOR SNIP_GRAD_FACES TOP */  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1=0; k1 < m_NE1; ++k1) {  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;  
                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;  
                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;  
                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;  
                     } /* end of component loop i */  
                 } /* end of k1 loop */  
             } /* end of face 0 */  
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1=0; k1 < m_NE1; ++k1) {  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;  
                         o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;  
                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;  
                         o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;  
                     } /* end of component loop i */  
                 } /* end of k1 loop */  
             } /* end of face 1 */  
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;  
                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;  
                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;  
                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of face 2 */  
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;  
                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;  
                         o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;  
                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of face 3 */  
             /* GENERATOR SNIP_GRAD_FACES BOTTOM */  
         } // end of parallel section  
     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
 #pragma omp parallel  
         {  
             /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1=0; k1 < m_NE1; ++k1) {  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);  
                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;  
                     } /* end of component loop i */  
                 } /* end of k1 loop */  
             } /* end of face 0 */  
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1=0; k1 < m_NE1; ++k1) {  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);  
                         o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;  
                     } /* end of component loop i */  
                 } /* end of k1 loop */  
             } /* end of face 1 */  
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;  
                         o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of face 2 */  
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0=0; k0 < m_NE0; ++k0) {  
                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));  
                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));  
                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));  
                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;  
                         o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);  
                     } /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of face 3 */  
             /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */  
         } // end of parallel section  
     } else {  
         stringstream msg;  
         msg << "setToGradient() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
 }  
   
 void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  
 {  
     escript::Data& in = *const_cast<escript::Data*>(&arg);  
     const dim_t numComp = in.getDataPointSize();  
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
     if (arg.getFunctionSpace().getTypeCode() == Elements) {  
         const double w_0 = h0*h1/4.;  
 #pragma omp parallel  
         {  
             vector<double> int_local(numComp, 0);  
 #pragma omp for nowait  
             for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         const register double f_2 = f[INDEX2(i,2,numComp)];  
                         const register double f_3 = f[INDEX2(i,3,numComp)];  
                         int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k1 loop */  
   
 #pragma omp critical  
             for (index_t i=0; i<numComp; i++)  
                 integrals[i]+=int_local[i];  
         } // end of parallel section  
     } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
         const double w_0 = h0*h1;  
 #pragma omp parallel  
         {  
             vector<double> int_local(numComp, 0);  
 #pragma omp for nowait  
             for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*w_0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             } /* end of k1 loop */  
   
 #pragma omp critical  
             for (index_t i=0; i<numComp; i++)  
                 integrals[i]+=int_local[i];  
         } // end of parallel section  
     } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
         const double w_0 = h0/2.;  
         const double w_1 = h1/2.;  
 #pragma omp parallel  
         {  
             vector<double> int_local(numComp, 0);  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         int_local[i]+=(f_0+f_1)*w_1;  
                     }  /* end of component loop i */  
                 } /* end of k1 loop */  
             }  
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         int_local[i]+=(f_0+f_1)*w_1;  
                     }  /* end of component loop i */  
                 } /* end of k1 loop */  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         int_local[i]+=(f_0+f_1)*w_0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         const register double f_0 = f[INDEX2(i,0,numComp)];  
                         const register double f_1 = f[INDEX2(i,1,numComp)];  
                         int_local[i]+=(f_0+f_1)*w_0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             }  
   
 #pragma omp critical  
             for (index_t i=0; i<numComp; i++)  
                 integrals[i]+=int_local[i];  
         } // end of parallel section  
     } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
 #pragma omp parallel  
         {  
             vector<double> int_local(numComp, 0);  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*h1;  
                     }  /* end of component loop i */  
                 } /* end of k1 loop */  
             }  
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*h1;  
                     }  /* end of component loop i */  
                 } /* end of k1 loop */  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*h0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);  
                     for (index_t i=0; i < numComp; ++i) {  
                         int_local[i]+=f[i]*h0;  
                     }  /* end of component loop i */  
                 } /* end of k0 loop */  
438              }              }
439            default:
440  #pragma omp critical              break;
             for (index_t i=0; i<numComp; i++)  
                 integrals[i]+=int_local[i];  
         } // end of parallel section  
     } else {  
         stringstream msg;  
         msg << "setToIntegrals() not implemented for "  
             << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
441      }      }
442    
443        stringstream msg;
444        msg << "ownSample: invalid function space type " << fsType;
445        throw RipleyException(msg.str());
446  }  }
447    
448  void Rectangle::setToNormal(escript::Data& out) const  void Rectangle::setToNormal(escript::Data& out) const
449  {  {
450      if (out.getFunctionSpace().getTypeCode() == FaceElements) {      if (out.getFunctionSpace().getTypeCode() == FaceElements) {
451            out.requireWrite();
452  #pragma omp parallel  #pragma omp parallel
453          {          {
454              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
# Line 683  void Rectangle::setToNormal(escript::Dat Line 500  void Rectangle::setToNormal(escript::Dat
500              }              }
501          } // end of parallel section          } // end of parallel section
502      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
503            out.requireWrite();
504  #pragma omp parallel  #pragma omp parallel
505          {          {
506              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
# Line 724  void Rectangle::setToNormal(escript::Dat Line 542  void Rectangle::setToNormal(escript::Dat
542    
543      } else {      } else {
544          stringstream msg;          stringstream msg;
545          msg << "setToNormal() not implemented for "          msg << "setToNormal: invalid function space type "
546              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());              << out.getFunctionSpace().getTypeCode();
547          throw RipleyException(msg.str());          throw RipleyException(msg.str());
548      }      }
549  }  }
550    
551  void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,  void Rectangle::setToSize(escript::Data& out) const
         escript::Data& rhs, const escript::Data& A, const escript::Data& B,  
         const escript::Data& C, const escript::Data& D,  
         const escript::Data& X, const escript::Data& Y,  
         const escript::Data& d, const escript::Data& y,  
         const escript::Data& d_contact, const escript::Data& y_contact,  
         const escript::Data& d_dirac, const escript::Data& y_dirac) const  
552  {  {
553      throw RipleyException("addPDEToSystem() not implemented");      if (out.getFunctionSpace().getTypeCode() == Elements
554                || out.getFunctionSpace().getTypeCode() == ReducedElements) {
555            out.requireWrite();
556            const dim_t numQuad=out.getNumDataPointsPerSample();
557            const double hSize=getFirstCoordAndSpacing(0).second;
558            const double vSize=getFirstCoordAndSpacing(1).second;
559            const double size=sqrt(hSize*hSize+vSize*vSize);
560    #pragma omp parallel for
561            for (index_t k = 0; k < getNumElements(); ++k) {
562                double* o = out.getSampleDataRW(k);
563                fill(o, o+numQuad, size);
564            }
565        } else if (out.getFunctionSpace().getTypeCode() == FaceElements
566                || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
567            out.requireWrite();
568            const dim_t numQuad=out.getNumDataPointsPerSample();
569            const double hSize=getFirstCoordAndSpacing(0).second;
570            const double vSize=getFirstCoordAndSpacing(1).second;
571    #pragma omp parallel
572            {
573                if (m_faceOffset[0] > -1) {
574    #pragma omp for nowait
575                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
576                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
577                        fill(o, o+numQuad, vSize);
578                    }
579                }
580    
581                if (m_faceOffset[1] > -1) {
582    #pragma omp for nowait
583                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
584                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
585                        fill(o, o+numQuad, vSize);
586                    }
587                }
588    
589                if (m_faceOffset[2] > -1) {
590    #pragma omp for nowait
591                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
592                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
593                        fill(o, o+numQuad, hSize);
594                    }
595                }
596    
597                if (m_faceOffset[3] > -1) {
598    #pragma omp for nowait
599                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
600                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
601                        fill(o, o+numQuad, hSize);
602                    }
603                }
604            } // end of parallel section
605    
606        } else {
607            stringstream msg;
608            msg << "setToSize: invalid function space type "
609                << out.getFunctionSpace().getTypeCode();
610            throw RipleyException(msg.str());
611        }
612  }  }
613    
614  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
615                                                  bool reducedColOrder) const                                                  bool reducedColOrder) const
616  {  {
617        /* FIXME: reduced
618      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
619          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
   
     // connector  
     RankVector neighbour;  
     IndexVector offsetInShared(1,0);  
     IndexVector sendShared, recvShared;  
     const IndexVector faces=getNumFacesPerBoundary();  
     const index_t left = (m_offset0==0 ? 0 : 1);  
     const index_t bottom = (m_offset1==0 ? 0 : 1);  
     // corner node from bottom-left  
     if (faces[0] == 0 && faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX-1);  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0+left]);  
         recvShared.push_back(m_nodeId[0]);  
     }  
     // bottom edge  
     if (faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX);  
         offsetInShared.push_back(offsetInShared.back()+m_N0-left);  
         for (dim_t i=left; i<m_N0; i++) {  
             // easy case, we know the neighbour id's  
             sendShared.push_back(m_nodeId[m_N0+i]);  
             recvShared.push_back(m_nodeId[i]);  
         }  
     }  
     // corner node from bottom-right  
     if (faces[1] == 0 && faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX+1);  
         const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour.back()/m_NX == 0 ? m_N1 : m_N1-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[(bottom+1)*m_N0-1]);  
         recvShared.push_back(first+N0*(N1-1));  
     }  
     // left edge  
     if (faces[0] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-1);  
         offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);  
         for (dim_t i=bottom; i<m_N1; i++) {  
             // easy case, we know the neighbour id's  
             sendShared.push_back(m_nodeId[i*m_N0+1]);  
             recvShared.push_back(m_nodeId[i*m_N0]);  
         }  
     }  
     // right edge  
     if (faces[1] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+1);  
         const index_t rightN0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);  
         for (dim_t i=bottom; i<m_N1; i++) {  
             sendShared.push_back(m_nodeId[(i+1)*m_N0-1]);  
             recvShared.push_back(first+rightN0*(i-bottom));  
         }  
     }  
     // corner node from top-left  
     if (faces[0] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX-1);  
         const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+left]);  
         recvShared.push_back(first+N0-1);  
     }  
     // top edge  
     if (faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+m_N0-left);  
         for (dim_t i=left; i<m_N0; i++) {  
             sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+i]);  
             recvShared.push_back(first+i-left);  
         }  
     }  
     // corner node from top-right  
     if (faces[1] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX+1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0*m_N1-1]);  
         recvShared.push_back(first);  
     }  
     const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];  
     /*  
     cout << "--- rcv_shcomp ---" << endl;  
     cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;  
     for (size_t i=0; i<neighbour.size(); i++) {  
         cout << "neighbor[" << i << "]=" << neighbour[i]  
             << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;  
     }  
     for (size_t i=0; i<recvShared.size(); i++) {  
         cout << "shared[" << i << "]=" << recvShared[i] << endl;  
     }  
     cout << "--- snd_shcomp ---" << endl;  
     for (size_t i=0; i<sendShared.size(); i++) {  
         cout << "shared[" << i << "]=" << sendShared[i] << endl;  
     }  
620      */      */
621        return m_pattern;
     Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &sendShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &recvShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);  
     Paso_SharedComponents_free(snd_shcomp);  
     Paso_SharedComponents_free(rcv_shcomp);  
   
     // create patterns  
     dim_t M, N;  
     IndexVector ptr(1,0);  
     IndexVector index;  
   
     // main pattern  
     for (index_t i=0; i<numDOF; i++) {  
         // always add the node itself  
         index.push_back(i);  
         int num=insertNeighbours(index, i);  
         ptr.push_back(ptr.back()+num+1);  
     }  
     M=N=ptr.size()-1;  
     // paso will manage the memory  
     index_t* indexC = MEMALLOC(index.size(),index_t);  
     index_t* ptrC = MEMALLOC(ptr.size(), index_t);  
     copy(index.begin(), index.end(), indexC);  
     copy(ptr.begin(), ptr.end(), ptrC);  
     Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
             M, N, ptrC, indexC);  
   
     /*  
     cout << "--- main_pattern ---" << endl;  
     cout << "M=" << M << ", N=" << N << endl;  
     for (size_t i=0; i<ptr.size(); i++) {  
         cout << "ptr[" << i << "]=" << ptr[i] << endl;  
     }  
     for (size_t i=0; i<index.size(); i++) {  
         cout << "index[" << i << "]=" << index[i] << endl;  
     }  
     */  
   
     ptr.clear();  
     index.clear();  
   
     // column & row couple patterns  
     Paso_Pattern *colCouplePattern, *rowCouplePattern;  
     generateCouplePatterns(&colCouplePattern, &rowCouplePattern);  
   
     // allocate paso distribution  
     Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,  
             const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);  
   
     Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(  
             MATRIX_FORMAT_DEFAULT, distribution, distribution,  
             mainPattern, colCouplePattern, rowCouplePattern,  
             connector, connector);  
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colCouplePattern);  
     Paso_Pattern_free(rowCouplePattern);  
     Paso_Distribution_free(distribution);  
     return pattern;  
622  }  }
623    
624  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
# Line 961  IndexVector Rectangle::getNumFacesPerBou Line 672  IndexVector Rectangle::getNumFacesPerBou
672      return ret;      return ret;
673  }  }
674    
675    IndexVector Rectangle::getNumSubdivisionsPerDim() const
676    {
677        IndexVector ret;
678        ret.push_back(m_NX);
679        ret.push_back(m_NY);
680        return ret;
681    }
682    
683  pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
684  {  {
685      if (dim==0) {      if (dim==0) {
686          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);
687      } else if (dim==1) {      } else if (dim==1) {
688          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);
689      }      }
690      throw RipleyException("getFirstCoordAndSpacing(): invalid argument");      throw RipleyException("getFirstCoordAndSpacing: invalid argument");
691  }  }
692    
693  //protected  //protected
694  dim_t Rectangle::getNumDOF() const  dim_t Rectangle::getNumDOF() const
695  {  {
696      return m_nodeDistribution[m_mpiInfo->rank+1]      return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;
         -m_nodeDistribution[m_mpiInfo->rank];  
697  }  }
698    
699  //protected  //protected
# Line 1011  void Rectangle::assembleCoordinates(escr Line 729  void Rectangle::assembleCoordinates(escr
729      }      }
730  }  }
731    
732  //private  //protected
733  void Rectangle::populateSampleIds()  void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
734  {  {
735      // identifiers are ordered from left to right, bottom to top on each rank,      const dim_t numComp = in.getDataPointSize();
736      // except for the shared nodes which are owned by the rank below / to the      const double h0 = m_l0/m_gNE0;
737      // left of the current rank      const double h1 = m_l1/m_gNE1;
738        const double cx0 = -1./h0;
739        const double cx1 = -.78867513459481288225/h0;
740        const double cx2 = -.5/h0;
741        const double cx3 = -.21132486540518711775/h0;
742        const double cx4 = .21132486540518711775/h0;
743        const double cx5 = .5/h0;
744        const double cx6 = .78867513459481288225/h0;
745        const double cx7 = 1./h0;
746        const double cy0 = -1./h1;
747        const double cy1 = -.78867513459481288225/h1;
748        const double cy2 = -.5/h1;
749        const double cy3 = -.21132486540518711775/h1;
750        const double cy4 = .21132486540518711775/h1;
751        const double cy5 = .5/h1;
752        const double cy6 = .78867513459481288225/h1;
753        const double cy7 = 1./h1;
754    
755      // build node distribution vector first.      if (out.getFunctionSpace().getTypeCode() == Elements) {
756      // m_nodeDistribution[i] is the first node id on rank i, that is          out.requireWrite();
757      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes  #pragma omp parallel
758      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);          {
759      m_nodeDistribution[1]=getNumNodes();              vector<double> f_00(numComp);
760      for (dim_t k=1; k<m_mpiInfo->size-1; k++) {              vector<double> f_01(numComp);
761          const index_t x=k%m_NX;              vector<double> f_10(numComp);
762          const index_t y=k/m_NX;              vector<double> f_11(numComp);
763          index_t numNodes=getNumNodes();  #pragma omp for
764          if (x>0)              for (index_t k1=0; k1 < m_NE1; ++k1) {
765              numNodes-=m_N1;                  for (index_t k0=0; k0 < m_NE0; ++k0) {
766          if (y>0)                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
767              numNodes-=m_N0;                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
768          if (x>0 && y>0)                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
769              numNodes++; // subtracted corner twice -> fix that                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
770          m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
771                        for (index_t i=0; i < numComp; ++i) {
772                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
773                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
774                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
775                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
776                            o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
777                            o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
778                            o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
779                            o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
780                        } // end of component loop i
781                    } // end of k0 loop
782                } // end of k1 loop
783            } // end of parallel section
784        } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
785            out.requireWrite();
786    #pragma omp parallel
787            {
788                vector<double> f_00(numComp);
789                vector<double> f_01(numComp);
790                vector<double> f_10(numComp);
791                vector<double> f_11(numComp);
792    #pragma omp for
793                for (index_t k1=0; k1 < m_NE1; ++k1) {
794                    for (index_t k0=0; k0 < m_NE0; ++k0) {
795                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
796                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
797                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
798                        memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
799                        double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
800                        for (index_t i=0; i < numComp; ++i) {
801                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
802                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
803                        } // end of component loop i
804                    } // end of k0 loop
805                } // end of k1 loop
806            } // end of parallel section
807        } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
808            out.requireWrite();
809    #pragma omp parallel
810            {
811                vector<double> f_00(numComp);
812                vector<double> f_01(numComp);
813                vector<double> f_10(numComp);
814                vector<double> f_11(numComp);
815                if (m_faceOffset[0] > -1) {
816    #pragma omp for nowait
817                    for (index_t k1=0; k1 < m_NE1; ++k1) {
818                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
819                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
820                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
821                        memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
822                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
823                        for (index_t i=0; i < numComp; ++i) {
824                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
825                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
826                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
827                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
828                        } // end of component loop i
829                    } // end of k1 loop
830                } // end of face 0
831                if (m_faceOffset[1] > -1) {
832    #pragma omp for nowait
833                    for (index_t k1=0; k1 < m_NE1; ++k1) {
834                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
835                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
836                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
837                        memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
838                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
839                        for (index_t i=0; i < numComp; ++i) {
840                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
841                            o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
842                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
843                            o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
844                        } // end of component loop i
845                    } // end of k1 loop
846                } // end of face 1
847                if (m_faceOffset[2] > -1) {
848    #pragma omp for nowait
849                    for (index_t k0=0; k0 < m_NE0; ++k0) {
850                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
851                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
852                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
853                        memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
854                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
855                        for (index_t i=0; i < numComp; ++i) {
856                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
857                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
858                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
859                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
860                        } // end of component loop i
861                    } // end of k0 loop
862                } // end of face 2
863                if (m_faceOffset[3] > -1) {
864    #pragma omp for nowait
865                    for (index_t k0=0; k0 < m_NE0; ++k0) {
866                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
867                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
868                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
869                        memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
870                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
871                        for (index_t i=0; i < numComp; ++i) {
872                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
873                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
874                            o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
875                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
876                        } // end of component loop i
877                    } // end of k0 loop
878                } // end of face 3
879            } // end of parallel section
880    
881        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
882            out.requireWrite();
883    #pragma omp parallel
884            {
885                vector<double> f_00(numComp);
886                vector<double> f_01(numComp);
887                vector<double> f_10(numComp);
888                vector<double> f_11(numComp);
889                if (m_faceOffset[0] > -1) {
890    #pragma omp for nowait
891                    for (index_t k1=0; k1 < m_NE1; ++k1) {
892                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
893                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
894                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
895                        memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
896                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
897                        for (index_t i=0; i < numComp; ++i) {
898                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
899                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
900                        } // end of component loop i
901                    } // end of k1 loop
902                } // end of face 0
903                if (m_faceOffset[1] > -1) {
904    #pragma omp for nowait
905                    for (index_t k1=0; k1 < m_NE1; ++k1) {
906                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
907                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
908                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
909                        memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
910                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
911                        for (index_t i=0; i < numComp; ++i) {
912                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
913                            o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
914                        } // end of component loop i
915                    } // end of k1 loop
916                } // end of face 1
917                if (m_faceOffset[2] > -1) {
918    #pragma omp for nowait
919                    for (index_t k0=0; k0 < m_NE0; ++k0) {
920                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
921                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
922                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
923                        memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
924                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
925                        for (index_t i=0; i < numComp; ++i) {
926                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
927                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
928                        } // end of component loop i
929                    } // end of k0 loop
930                } // end of face 2
931                if (m_faceOffset[3] > -1) {
932    #pragma omp for nowait
933                    for (index_t k0=0; k0 < m_NE0; ++k0) {
934                        memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
935                        memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
936                        memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
937                        memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
938                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
939                        for (index_t i=0; i < numComp; ++i) {
940                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
941                            o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
942                        } // end of component loop i
943                    } // end of k0 loop
944                } // end of face 3
945            } // end of parallel section
946      }      }
947      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();  }
948    
949      m_dofId.resize(getNumDOF());  //protected
950      m_nodeId.resize(getNumNodes());  void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
951    {
952        const dim_t numComp = arg.getDataPointSize();
953        const double h0 = m_l0/m_gNE0;
954        const double h1 = m_l1/m_gNE1;
955        const index_t left = (m_offset0==0 ? 0 : 1);
956        const index_t bottom = (m_offset1==0 ? 0 : 1);
957        const int fs=arg.getFunctionSpace().getTypeCode();
958        if (fs == Elements && arg.actsExpanded()) {
959    #pragma omp parallel
960            {
961                vector<double> int_local(numComp, 0);
962                const double w = h0*h1/4.;
963    #pragma omp for nowait
964                for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
965                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
966                        const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
967                        for (index_t i=0; i < numComp; ++i) {
968                            const double f0 = f[INDEX2(i,0,numComp)];
969                            const double f1 = f[INDEX2(i,1,numComp)];
970                            const double f2 = f[INDEX2(i,2,numComp)];
971                            const double f3 = f[INDEX2(i,3,numComp)];
972                            int_local[i]+=(f0+f1+f2+f3)*w;
973                        }  // end of component loop i
974                    } // end of k0 loop
975                } // end of k1 loop
976    #pragma omp critical
977                for (index_t i=0; i<numComp; i++)
978                    integrals[i]+=int_local[i];
979            } // end of parallel section
980    
981        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
982            const double w = h0*h1;
983    #pragma omp parallel
984            {
985                vector<double> int_local(numComp, 0);
986    #pragma omp for nowait
987                for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
988                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
989                        const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
990                        for (index_t i=0; i < numComp; ++i) {
991                            int_local[i]+=f[i]*w;
992                        }
993                    }
994                }
995    #pragma omp critical
996                for (index_t i=0; i<numComp; i++)
997                    integrals[i]+=int_local[i];
998            } // end of parallel section
999    
1000        } else if (fs == FaceElements && arg.actsExpanded()) {
1001    #pragma omp parallel
1002            {
1003                vector<double> int_local(numComp, 0);
1004                const double w0 = h0/2.;
1005                const double w1 = h1/2.;
1006                if (m_faceOffset[0] > -1) {
1007    #pragma omp for nowait
1008                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1009                        const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1010                        for (index_t i=0; i < numComp; ++i) {
1011                            const double f0 = f[INDEX2(i,0,numComp)];
1012                            const double f1 = f[INDEX2(i,1,numComp)];
1013                            int_local[i]+=(f0+f1)*w1;
1014                        }  // end of component loop i
1015                    } // end of k1 loop
1016                }
1017    
1018                if (m_faceOffset[1] > -1) {
1019    #pragma omp for nowait
1020                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1021                        const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1022                        for (index_t i=0; i < numComp; ++i) {
1023                            const double f0 = f[INDEX2(i,0,numComp)];
1024                            const double f1 = f[INDEX2(i,1,numComp)];
1025                            int_local[i]+=(f0+f1)*w1;
1026                        }  // end of component loop i
1027                    } // end of k1 loop
1028                }
1029    
1030                if (m_faceOffset[2] > -1) {
1031    #pragma omp for nowait
1032                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1033                        const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1034                        for (index_t i=0; i < numComp; ++i) {
1035                            const double f0 = f[INDEX2(i,0,numComp)];
1036                            const double f1 = f[INDEX2(i,1,numComp)];
1037                            int_local[i]+=(f0+f1)*w0;
1038                        }  // end of component loop i
1039                    } // end of k0 loop
1040                }
1041    
1042                if (m_faceOffset[3] > -1) {
1043    #pragma omp for nowait
1044                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1045                        const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1046                        for (index_t i=0; i < numComp; ++i) {
1047                            const double f0 = f[INDEX2(i,0,numComp)];
1048                            const double f1 = f[INDEX2(i,1,numComp)];
1049                            int_local[i]+=(f0+f1)*w0;
1050                        }  // end of component loop i
1051                    } // end of k0 loop
1052                }
1053    #pragma omp critical
1054                for (index_t i=0; i<numComp; i++)
1055                    integrals[i]+=int_local[i];
1056            } // end of parallel section
1057    
1058        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1059    #pragma omp parallel
1060            {
1061                vector<double> int_local(numComp, 0);
1062                if (m_faceOffset[0] > -1) {
1063    #pragma omp for nowait
1064                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1065                        const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1066                        for (index_t i=0; i < numComp; ++i) {
1067                            int_local[i]+=f[i]*h1;
1068                        }
1069                    }
1070                }
1071    
1072                if (m_faceOffset[1] > -1) {
1073    #pragma omp for nowait
1074                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1075                        const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1076                        for (index_t i=0; i < numComp; ++i) {
1077                            int_local[i]+=f[i]*h1;
1078                        }
1079                    }
1080                }
1081    
1082                if (m_faceOffset[2] > -1) {
1083    #pragma omp for nowait
1084                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1085                        const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1086                        for (index_t i=0; i < numComp; ++i) {
1087                            int_local[i]+=f[i]*h0;
1088                        }
1089                    }
1090                }
1091    
1092                if (m_faceOffset[3] > -1) {
1093    #pragma omp for nowait
1094                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1095                        const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1096                        for (index_t i=0; i < numComp; ++i) {
1097                            int_local[i]+=f[i]*h0;
1098                        }
1099                    }
1100                }
1101    
1102    #pragma omp critical
1103                for (index_t i=0; i<numComp; i++)
1104                    integrals[i]+=int_local[i];
1105            } // end of parallel section
1106        } // function space selector
1107    }
1108    
1109    //protected
1110    dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
1111    {
1112        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1113        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1114        const int x=node%nDOF0;
1115        const int y=node/nDOF0;
1116        dim_t num=0;
1117        // loop through potential neighbours and add to index if positions are
1118        // within bounds
1119        for (int i1=-1; i1<2; i1++) {
1120            for (int i0=-1; i0<2; i0++) {
1121                // skip node itself
1122                if (i0==0 && i1==0)
1123                    continue;
1124                // location of neighbour node
1125                const int nx=x+i0;
1126                const int ny=y+i1;
1127                if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {
1128                    index.push_back(ny*nDOF0+nx);
1129                    num++;
1130                }
1131            }
1132        }
1133    
1134        return num;
1135    }
1136    
1137    //protected
1138    void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const
1139    {
1140        const dim_t numComp = in.getDataPointSize();
1141        out.requireWrite();
1142    
     // the bottom row and left column are not owned by this rank so the  
     // identifiers need to be computed accordingly  
1143      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1144      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1145      if (left>0) {      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1146          const int neighbour=m_mpiInfo->rank-1;      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
         const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
1147  #pragma omp parallel for  #pragma omp parallel for
1148          for (dim_t i1=bottom; i1<m_N1; i1++) {      for (index_t i=0; i<nDOF1; i++) {
1149              m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]          for (index_t j=0; j<nDOF0; j++) {
1150                  + (i1-bottom+1)*leftN0-1;              const index_t n=j+left+(i+bottom)*m_N0;
1151                const double* src=in.getSampleDataRO(n);
1152                copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
1153          }          }
1154      }      }
1155      if (bottom>0) {  }
1156          const int neighbour=m_mpiInfo->rank-m_NX;  
1157          const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  //protected
1158          const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);  void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const
1159    {
1160        const dim_t numComp = in.getDataPointSize();
1161        Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
1162        in.requireWrite();
1163        Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));
1164    
1165        const dim_t numDOF = getNumDOF();
1166        out.requireWrite();
1167        const double* buffer = Paso_Coupler_finishCollect(coupler);
1168    
1169  #pragma omp parallel for  #pragma omp parallel for
1170          for (dim_t i0=left; i0<m_N0; i0++) {      for (index_t i=0; i<getNumNodes(); i++) {
1171              m_nodeId[i0]=m_nodeDistribution[neighbour]          const double* src=(m_dofMap[i]<numDOF ?
1172                  + (bottomN1-1)*bottomN0 + i0 - left;                  in.getSampleDataRO(m_dofMap[i])
1173          }                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1174            copy(src, src+numComp, out.getSampleDataRW(i));
1175      }      }
1176      if (left>0 && bottom>0) {      Paso_Coupler_free(coupler);
1177          const int neighbour=m_mpiInfo->rank-m_NX-1;  }
1178          const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
1179          const index_t N1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);  //private
1180          m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1-1;  void Rectangle::populateSampleIds()
1181    {
1182        // identifiers are ordered from left to right, bottom to top globablly.
1183    
1184        // build node distribution vector first.
1185        // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
1186        m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1187        const dim_t numDOF=getNumDOF();
1188        for (dim_t k=1; k<m_mpiInfo->size; k++) {
1189            m_nodeDistribution[k]=k*numDOF;
1190      }      }
1191        m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
1192        m_nodeId.resize(getNumNodes());
1193        m_dofId.resize(numDOF);
1194        m_elementId.resize(getNumElements());
1195        m_faceId.resize(getNumFaceElements());
1196    
1197      // the rest of the id's are contiguous  #pragma omp parallel
1198      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];      {
1199  #pragma omp parallel for          // nodes
1200      for (dim_t i1=bottom; i1<m_N1; i1++) {  #pragma omp for nowait
1201          for (dim_t i0=left; i0<m_N0; i0++) {          for (dim_t i1=0; i1<m_N1; i1++) {
1202              const index_t idx=i0-left+(i1-bottom)*(m_N0-left);              for (dim_t i0=0; i0<m_N0; i0++) {
1203              m_nodeId[i0+i1*m_N0] = firstId+idx;                  m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;
1204              m_dofId[idx] = firstId+idx;              }
1205          }          }
1206      }  
1207            // degrees of freedom
1208    #pragma omp for nowait
1209            for (dim_t k=0; k<numDOF; k++)
1210                m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;
1211    
1212            // elements
1213    #pragma omp for nowait
1214            for (dim_t i1=0; i1<m_NE1; i1++) {
1215                for (dim_t i0=0; i0<m_NE0; i0++) {
1216                    m_elementId[i0+i1*m_NE0]=(m_offset1+i1)*m_gNE0+m_offset0+i0;
1217                }
1218            }
1219    
1220            // face elements
1221    #pragma omp for
1222            for (dim_t k=0; k<getNumFaceElements(); k++)
1223                m_faceId[k]=k;
1224        } // end parallel section
1225    
1226      m_nodeTags.assign(getNumNodes(), 0);      m_nodeTags.assign(getNumNodes(), 0);
1227      updateTagsInUse(Nodes);      updateTagsInUse(Nodes);
1228    
     // elements  
     m_elementId.resize(getNumElements());  
 #pragma omp parallel for  
     for (dim_t k=0; k<getNumElements(); k++) {  
         m_elementId[k]=k;  
     }  
1229      m_elementTags.assign(getNumElements(), 0);      m_elementTags.assign(getNumElements(), 0);
1230      updateTagsInUse(Elements);      updateTagsInUse(Elements);
1231    
     // face element id's  
     m_faceId.resize(getNumFaceElements());  
 #pragma omp parallel for  
     for (dim_t k=0; k<getNumFaceElements(); k++) {  
         m_faceId[k]=k;  
     }  
   
1232      // generate face offset vector and set face tags      // generate face offset vector and set face tags
1233      const IndexVector facesPerEdge = getNumFacesPerBoundary();      const IndexVector facesPerEdge = getNumFacesPerBoundary();
1234      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
# Line 1121  void Rectangle::populateSampleIds() Line 1251  void Rectangle::populateSampleIds()
1251  }  }
1252    
1253  //private  //private
1254  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const  void Rectangle::createPattern()
1255  {  {
1256      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1257      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1258      const int x=node%myN0;      const index_t left = (m_offset0==0 ? 0 : 1);
1259      const int y=node/myN0;      const index_t bottom = (m_offset1==0 ? 0 : 1);
1260      int num=0;  
1261      if (y>0) {      // populate node->DOF mapping with own degrees of freedom.
1262          if (x>0) {      // The rest is assigned in the loop further down
1263              // bottom-left      m_dofMap.assign(getNumNodes(), 0);
1264              index.push_back(node-myN0-1);  #pragma omp parallel for
1265              num++;      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1266          }          for (index_t j=left; j<left+nDOF0; j++) {
1267          // bottom              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
1268          index.push_back(node-myN0);          }
         num++;  
         if (x<myN0-1) {  
             // bottom-right  
             index.push_back(node-myN0+1);  
             num++;  
         }  
     }  
     if (x<myN0-1) {  
         // right  
         index.push_back(node+1);  
         num++;  
         if (y<myN1-1) {  
             // top-right  
             index.push_back(node+myN0+1);  
             num++;  
         }  
     }  
     if (y<myN1-1) {  
         // top  
         index.push_back(node+myN0);  
         num++;  
         if (x>0) {  
             // top-left  
             index.push_back(node+myN0-1);  
             num++;  
         }  
     }  
     if (x>0) {  
         // left  
         index.push_back(node-1);  
         num++;  
1269      }      }
1270    
1271      return num;      // build list of shared components and neighbours by looping through
1272  }      // all potential neighbouring ranks and checking if positions are
1273        // within bounds
1274        const dim_t numDOF=nDOF0*nDOF1;
1275        vector<IndexVector> colIndices(numDOF); // for the couple blocks
1276        RankVector neighbour;
1277        IndexVector offsetInShared(1,0);
1278        IndexVector sendShared, recvShared;
1279        int numShared=0;
1280        const int x=m_mpiInfo->rank%m_NX;
1281        const int y=m_mpiInfo->rank/m_NX;
1282        for (int i1=-1; i1<2; i1++) {
1283            for (int i0=-1; i0<2; i0++) {
1284                // skip this rank
1285                if (i0==0 && i1==0)
1286                    continue;
1287                // location of neighbour rank
1288                const int nx=x+i0;
1289                const int ny=y+i1;
1290                if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {
1291                    neighbour.push_back(ny*m_NX+nx);
1292                    if (i0==0) {
1293                        // sharing top or bottom edge
1294                        const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1295                        const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);
1296                        offsetInShared.push_back(offsetInShared.back()+nDOF0);
1297                        for (dim_t i=0; i<nDOF0; i++, numShared++) {
1298                            sendShared.push_back(firstDOF+i);
1299                            recvShared.push_back(numDOF+numShared);
1300                            if (i>0)
1301                                colIndices[firstDOF+i-1].push_back(numShared);
1302                            colIndices[firstDOF+i].push_back(numShared);
1303                            if (i<nDOF0-1)
1304                                colIndices[firstDOF+i+1].push_back(numShared);
1305                            m_dofMap[firstNode+i]=numDOF+numShared;
1306                        }
1307                    } else if (i1==0) {
1308                        // sharing left or right edge
1309                        const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1310                        const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);
1311                        offsetInShared.push_back(offsetInShared.back()+nDOF1);
1312                        for (dim_t i=0; i<nDOF1; i++, numShared++) {
1313                            sendShared.push_back(firstDOF+i*nDOF0);
1314                            recvShared.push_back(numDOF+numShared);
1315                            if (i>0)
1316                                colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);
1317                            colIndices[firstDOF+i*nDOF0].push_back(numShared);
1318                            if (i<nDOF1-1)
1319                                colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1320                            m_dofMap[firstNode+i*m_N0]=numDOF+numShared;
1321                        }
1322                    } else {
1323                        // sharing a node
1324                        const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1325                        const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);
1326                        offsetInShared.push_back(offsetInShared.back()+1);
1327                        sendShared.push_back(dof);
1328                        recvShared.push_back(numDOF+numShared);
1329                        colIndices[dof].push_back(numShared);
1330                        m_dofMap[node]=numDOF+numShared;
1331                        ++numShared;
1332                    }
1333                }
1334            }
1335        }
1336    
1337  //private      // create connector
1338  void Rectangle::generateCouplePatterns(Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1339  {              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
1340      IndexVector ptr(1,0);              &offsetInShared[0], 1, 0, m_mpiInfo);
1341      IndexVector index;      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
1342      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
1343      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);              &offsetInShared[0], 1, 0, m_mpiInfo);
1344      const IndexVector faces=getNumFacesPerBoundary();      m_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1345        Paso_SharedComponents_free(snd_shcomp);
1346        Paso_SharedComponents_free(rcv_shcomp);
1347    
1348      // bottom edge      // create main and couple blocks
1349      dim_t n=0;      Paso_Pattern *mainPattern = createMainPattern();
1350      if (faces[0] == 0) {      Paso_Pattern *colPattern, *rowPattern;
1351          index.push_back(2*(myN0+myN1+1));      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
1352          index.push_back(2*(myN0+myN1+1)+1);  
1353          n+=2;      // allocate paso distribution
1354          if (faces[2] == 0) {      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
1355              index.push_back(0);              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
1356              index.push_back(1);  
1357              index.push_back(2);      // finally create the system matrix
1358              n+=3;      m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
1359          }              distribution, distribution, mainPattern, colPattern, rowPattern,
1360      } else if (faces[2] == 0) {              m_connector, m_connector);
1361          index.push_back(1);  
1362          index.push_back(2);      Paso_Distribution_free(distribution);
1363          n+=2;  
1364      }      // useful debug output
1365      // n=neighbours of bottom-left corner node      /*
1366      ptr.push_back(ptr.back()+n);      cout << "--- rcv_shcomp ---" << endl;
1367      n=0;      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
1368      if (faces[2] == 0) {      for (size_t i=0; i<neighbour.size(); i++) {
1369          for (dim_t i=1; i<myN0-1; i++) {          cout << "neighbor[" << i << "]=" << neighbour[i]
1370              index.push_back(i);              << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
1371              index.push_back(i+1);      }
1372              index.push_back(i+2);      for (size_t i=0; i<recvShared.size(); i++) {
1373              ptr.push_back(ptr.back()+3);          cout << "shared[" << i << "]=" << recvShared[i] << endl;
         }  
         index.push_back(myN0-1);  
         index.push_back(myN0);  
         n+=2;  
         if (faces[1] == 0) {  
             index.push_back(myN0+1);  
             index.push_back(myN0+2);  
             index.push_back(myN0+3);  
             n+=3;  
         }  
     } else {  
         for (dim_t i=1; i<myN0-1; i++) {  
             ptr.push_back(ptr.back());  
         }  
         if (faces[1] == 0) {  
             index.push_back(myN0+2);  
             index.push_back(myN0+3);  
             n+=2;  
         }  
     }  
     // n=neighbours of bottom-right corner node  
     ptr.push_back(ptr.back()+n);  
   
     // 2nd row to 2nd last row  
     for (dim_t i=1; i<myN1-1; i++) {  
         // left edge  
         if (faces[0] == 0) {  
             index.push_back(2*(myN0+myN1+2)-i);  
             index.push_back(2*(myN0+myN1+2)-i-1);  
             index.push_back(2*(myN0+myN1+2)-i-2);  
             ptr.push_back(ptr.back()+3);  
         } else {  
             ptr.push_back(ptr.back());  
         }  
         for (dim_t j=1; j<myN0-1; j++) {  
             ptr.push_back(ptr.back());  
         }  
         // right edge  
         if (faces[1] == 0) {  
             index.push_back(myN0+i+1);  
             index.push_back(myN0+i+2);  
             index.push_back(myN0+i+3);  
             ptr.push_back(ptr.back()+3);  
         } else {  
             ptr.push_back(ptr.back());  
         }  
     }  
   
     // top edge  
     n=0;  
     if (faces[0] == 0) {  
         index.push_back(2*myN0+myN1+5);  
         index.push_back(2*myN0+myN1+4);  
         n+=2;  
         if (faces[3] == 0) {  
             index.push_back(2*myN0+myN1+3);  
             index.push_back(2*myN0+myN1+2);  
             index.push_back(2*myN0+myN1+1);  
             n+=3;  
         }  
     } else if (faces[3] == 0) {  
         index.push_back(2*myN0+myN1+2);  
         index.push_back(2*myN0+myN1+1);  
         n+=2;  
     }  
     // n=neighbours of top-left corner node  
     ptr.push_back(ptr.back()+n);  
     n=0;  
     if (faces[3] == 0) {  
         for (dim_t i=1; i<myN0-1; i++) {  
             index.push_back(2*myN0+myN1+i+1);  
             index.push_back(2*myN0+myN1+i);  
             index.push_back(2*myN0+myN1+i-1);  
             ptr.push_back(ptr.back()+3);  
         }  
         index.push_back(myN0+myN1+4);  
         index.push_back(myN0+myN1+3);  
         n+=2;  
         if (faces[1] == 0) {  
             index.push_back(myN0+myN1+2);  
             index.push_back(myN0+myN1+1);  
             index.push_back(myN0+myN1);  
             n+=3;  
         }  
     } else {  
         for (dim_t i=1; i<myN0-1; i++) {  
             ptr.push_back(ptr.back());  
         }  
         if (faces[1] == 0) {  
             index.push_back(myN0+myN1+1);  
             index.push_back(myN0+myN1);  
             n+=2;  
         }  
     }  
     // n=neighbours of top-right corner node  
     ptr.push_back(ptr.back()+n);  
   
     dim_t M=ptr.size()-1;  
     map<index_t,index_t> idMap;  
     dim_t N=0;  
     for (IndexVector::iterator it=index.begin(); it!=index.end(); it++) {  
         if (idMap.find(*it)==idMap.end()) {  
             idMap[*it]=N;  
             *it=N++;  
         } else {  
             *it=idMap[*it];  
         }  
1374      }      }
1375        cout << "--- snd_shcomp ---" << endl;
1376        for (size_t i=0; i<sendShared.size(); i++) {
1377            cout << "shared[" << i << "]=" << sendShared[i] << endl;
1378        }
1379        cout << "--- dofMap ---" << endl;
1380        for (size_t i=0; i<m_dofMap.size(); i++) {
1381            cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
1382        }
1383        cout << "--- colIndices ---" << endl;
1384        for (size_t i=0; i<colIndices.size(); i++) {
1385            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
1386        }
1387        */
1388    
1389      /*      /*
1390      cout << "--- colCouple_pattern ---" << endl;      cout << "--- main_pattern ---" << endl;
1391      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
1392      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<mainPattern->numOutput+1; i++) {
1393          cout << "ptr[" << i << "]=" << ptr[i] << endl;          cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
1394      }      }
1395      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
1396          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
1397      }      }
1398      */      */
1399    
1400      // now build the row couple pattern      /*
1401      IndexVector ptr2(1,0);      cout << "--- colCouple_pattern ---" << endl;
1402      IndexVector index2;      cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
1403      for (dim_t id=0; id<N; id++) {      for (size_t i=0; i<colPattern->numOutput+1; i++) {
1404          n=0;          cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
         for (dim_t i=0; i<M; i++) {  
             for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {  
                 if (index[j]==id) {  
                     index2.push_back(i);  
                     n++;  
                     break;  
                 }  
             }  
         }  
         ptr2.push_back(ptr2.back()+n);  
1405      }      }
1406        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
1407            cout << "index[" << i << "]=" << colPattern->index[i] << endl;
1408        }
1409        */
1410    
1411      /*      /*
1412      cout << "--- rowCouple_pattern ---" << endl;      cout << "--- rowCouple_pattern ---" << endl;
1413      cout << "M=" << N << ", N=" << M << endl;      cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
1414      for (size_t i=0; i<ptr2.size(); i++) {      for (size_t i=0; i<rowPattern->numOutput+1; i++) {
1415          cout << "ptr[" << i << "]=" << ptr2[i] << endl;          cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
1416      }      }
1417      for (size_t i=0; i<index2.size(); i++) {      for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
1418          cout << "index[" << i << "]=" << index2[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
1419      }      }
1420      */      */
1421    
1422      // paso will manage the memory      Paso_Pattern_free(mainPattern);
1423      index_t* indexC = MEMALLOC(index.size(), index_t);      Paso_Pattern_free(colPattern);
1424      index_t* ptrC = MEMALLOC(ptr.size(), index_t);      Paso_Pattern_free(rowPattern);
1425      copy(index.begin(), index.end(), indexC);  }
1426      copy(ptr.begin(), ptr.end(), ptrC);  
1427      *colPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);  //private
1428    void Rectangle::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
1429      // paso will manage the memory           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
1430      indexC = MEMALLOC(index2.size(), index_t);           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
1431      ptrC = MEMALLOC(ptr2.size(), index_t);  {
1432      copy(index2.begin(), index2.end(), indexC);      IndexVector rowIndex;
1433      copy(ptr2.begin(), ptr2.end(), ptrC);      rowIndex.push_back(m_dofMap[firstNode]);
1434      *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);      rowIndex.push_back(m_dofMap[firstNode+1]);
1435        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
1436        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
1437        if (addF) {
1438            double *F_p=F.getSampleDataRW(0);
1439            for (index_t i=0; i<rowIndex.size(); i++) {
1440                if (rowIndex[i]<getNumDOF()) {
1441                    for (index_t eq=0; eq<nEq; eq++) {
1442                        F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
1443                    }
1444                }
1445            }
1446        }
1447        if (addS) {
1448            addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
1449        }
1450  }  }
1451    
1452  //protected  //protected
# Line 1376  void Rectangle::interpolateNodesOnElemen Line 1455  void Rectangle::interpolateNodesOnElemen
1455  {  {
1456      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1457      if (reduced) {      if (reduced) {
1458          /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */          out.requireWrite();
1459          const double c0 = .25;          const double c0 = 0.25;
1460  #pragma omp parallel for  #pragma omp parallel
1461          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1462              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1463                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_01(numComp);
1464                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1465                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              vector<double> f_11(numComp);
1466                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1467                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1468                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1469                      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));
1470                  } /* end of component loop i */                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
1471              } /* end of k0 loop */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
1472          } /* end of k1 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
1473          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1474                        for (index_t i=0; i < numComp; ++i) {
1475                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1476                        } /* end of component loop i */
1477                    } /* end of k0 loop */
1478                } /* end of k1 loop */
1479            } /* end of parallel section */
1480      } else {      } else {
1481          /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */          out.requireWrite();
1482          const double c0 = .16666666666666666667;          const double c0 = 0.16666666666666666667;
1483          const double c1 = .044658198738520451079;          const double c1 = 0.044658198738520451079;
1484          const double c2 = .62200846792814621559;          const double c2 = 0.62200846792814621559;
1485  #pragma omp parallel for  #pragma omp parallel
1486          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1487              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1488                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
1489                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1490                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
1491                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1492                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1493                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1494                      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));
1495                      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));
1496                      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));
1497                      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));
1498                  } /* end of component loop i */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1499              } /* end of k0 loop */                      for (index_t i=0; i < numComp; ++i) {
1500          } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1501          /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */                          o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1502                            o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1503                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1504                        } /* end of component loop i */
1505                    } /* end of k0 loop */
1506                } /* end of k1 loop */
1507            } /* end of parallel section */
1508      }      }
1509  }  }
1510    
# Line 1423  void Rectangle::interpolateNodesOnFaces( Line 1514  void Rectangle::interpolateNodesOnFaces(
1514  {  {
1515      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1516      if (reduced) {      if (reduced) {
1517          const double c0 = .5;          out.requireWrite();
1518            const double c0 = 0.5;
1519  #pragma omp parallel  #pragma omp parallel
1520          {          {
1521              /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */              vector<double> f_00(numComp);
1522                vector<double> f_01(numComp);
1523                vector<double> f_10(numComp);
1524                vector<double> f_11(numComp);
1525              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1526  #pragma omp for nowait  #pragma omp for nowait
1527                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1528                      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));
1529                      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));
1530                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1531                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1532                          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 1441  void Rectangle::interpolateNodesOnFaces( Line 1536  void Rectangle::interpolateNodesOnFaces(
1536              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1537  #pragma omp for nowait  #pragma omp for nowait
1538                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1539                      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));
1540                      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));
1541                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1542                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1543                          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 1452  void Rectangle::interpolateNodesOnFaces( Line 1547  void Rectangle::interpolateNodesOnFaces(
1547              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1548  #pragma omp for nowait  #pragma omp for nowait
1549                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1550                      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));
1551                      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));
1552                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1553                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1554                          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 1463  void Rectangle::interpolateNodesOnFaces( Line 1558  void Rectangle::interpolateNodesOnFaces(
1558              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1559  #pragma omp for nowait  #pragma omp for nowait
1560                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1561                      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));
1562                      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));
1563                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1564                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1565                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1566                      } /* end of component loop i */                      } /* end of component loop i */
1567                  } /* end of k0 loop */                  } /* end of k0 loop */
1568              } /* end of face 3 */              } /* end of face 3 */
1569              /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */          } /* end of parallel section */
         } // end of parallel section  
1570      } else {      } else {
1571            out.requireWrite();
1572          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1573          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1574  #pragma omp parallel  #pragma omp parallel
1575          {          {
1576              /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */              vector<double> f_00(numComp);
1577                vector<double> f_01(numComp);
1578                vector<double> f_10(numComp);
1579                vector<double> f_11(numComp);
1580              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1581  #pragma omp for nowait  #pragma omp for nowait
1582                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1583                      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));
1584                      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));
1585                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1586                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1587                          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];
1588                          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];
1589                      } /* end of component loop i */                      } /* end of component loop i */
1590                  } /* end of k1 loop */                  } /* end of k1 loop */
1591              } /* end of face 0 */              } /* end of face 0 */
1592              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1593  #pragma omp for nowait  #pragma omp for nowait
1594                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1595                      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));
1596                      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));
1597                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1598                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1599                          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];
1600                          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];
1601                      } /* end of component loop i */                      } /* end of component loop i */
1602                  } /* end of k1 loop */                  } /* end of k1 loop */
1603              } /* end of face 1 */              } /* end of face 1 */
1604              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1605  #pragma omp for nowait  #pragma omp for nowait
1606                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1607                      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));
1608                      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));
1609                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1610                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1611                          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];
1612                          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];
1613                      } /* end of component loop i */                      } /* end of component loop i */
1614                  } /* end of k0 loop */                  } /* end of k0 loop */
1615              } /* end of face 2 */              } /* end of face 2 */
1616              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1617  #pragma omp for nowait  #pragma omp for nowait
1618                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1619                      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));
1620                      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));
1621                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1622                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1623                          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];
1624                          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];
1625                      } /* end of component loop i */                      } /* end of component loop i */
1626                  } /* end of k0 loop */                  } /* end of k0 loop */
1627              } /* end of face 3 */              } /* end of face 3 */
1628              /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */          } /* end of parallel section */
1629          } // end of parallel section      }
1630    }
1631    
1632    //protected
1633    void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1634            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1635            const escript::Data& C, const escript::Data& D,
1636            const escript::Data& X, const escript::Data& Y) const
1637    {
1638        const double h0 = m_l0/m_gNE0;
1639        const double h1 = m_l1/m_gNE1;
1640        const double w0 = -0.1555021169820365539*h1/h0;
1641        const double w1 = 0.041666666666666666667;
1642        const double w2 = -0.15550211698203655390;
1643        const double w3 = 0.041666666666666666667*h0/h1;
1644        const double w4 = 0.15550211698203655390;
1645        const double w5 = -0.041666666666666666667;
1646        const double w6 = -0.01116454968463011277*h1/h0;
1647        const double w7 = 0.011164549684630112770;
1648        const double w8 = -0.011164549684630112770;
1649        const double w9 = -0.041666666666666666667*h1/h0;
1650        const double w10 = -0.041666666666666666667*h0/h1;
1651        const double w11 = 0.1555021169820365539*h1/h0;
1652        const double w12 = 0.1555021169820365539*h0/h1;
1653        const double w13 = 0.01116454968463011277*h0/h1;
1654        const double w14 = 0.01116454968463011277*h1/h0;
1655        const double w15 = 0.041666666666666666667*h1/h0;
1656        const double w16 = -0.01116454968463011277*h0/h1;
1657        const double w17 = -0.1555021169820365539*h0/h1;
1658        const double w18 = -0.33333333333333333333*h1/h0;
1659        const double w19 = 0.25;
1660        const double w20 = -0.25;
1661        const double w21 = 0.16666666666666666667*h0/h1;
1662        const double w22 = -0.16666666666666666667*h1/h0;
1663        const double w23 = -0.16666666666666666667*h0/h1;
1664        const double w24 = 0.33333333333333333333*h1/h0;
1665        const double w25 = 0.33333333333333333333*h0/h1;
1666        const double w26 = 0.16666666666666666667*h1/h0;
1667        const double w27 = -0.33333333333333333333*h0/h1;
1668        const double w28 = -0.032861463941450536761*h1;
1669        const double w29 = -0.032861463941450536761*h0;
1670        const double w30 = -0.12264065304058601714*h1;
1671        const double w31 = -0.0023593469594139828636*h1;
1672        const double w32 = -0.008805202725216129906*h0;
1673        const double w33 = -0.008805202725216129906*h1;
1674        const double w34 = 0.032861463941450536761*h1;
1675        const double w35 = 0.008805202725216129906*h1;
1676        const double w36 = 0.008805202725216129906*h0;
1677        const double w37 = 0.0023593469594139828636*h1;
1678        const double w38 = 0.12264065304058601714*h1;
1679        const double w39 = 0.032861463941450536761*h0;
1680        const double w40 = -0.12264065304058601714*h0;
1681        const double w41 = -0.0023593469594139828636*h0;
1682        const double w42 = 0.0023593469594139828636*h0;
1683        const double w43 = 0.12264065304058601714*h0;
1684        const double w44 = -0.16666666666666666667*h1;
1685        const double w45 = -0.083333333333333333333*h0;
1686        const double w46 = 0.083333333333333333333*h1;
1687        const double w47 = 0.16666666666666666667*h1;
1688        const double w48 = 0.083333333333333333333*h0;
1689        const double w49 = -0.16666666666666666667*h0;
1690        const double w50 = 0.16666666666666666667*h0;
1691        const double w51 = -0.083333333333333333333*h1;
1692        const double w52 = 0.025917019497006092316*h0*h1;
1693        const double w53 = 0.0018607582807716854616*h0*h1;
1694        const double w54 = 0.0069444444444444444444*h0*h1;
1695        const double w55 = 0.09672363354357992482*h0*h1;
1696        const double w56 = 0.00049858867864229740201*h0*h1;
1697        const double w57 = 0.055555555555555555556*h0*h1;
1698        const double w58 = 0.027777777777777777778*h0*h1;
1699        const double w59 = 0.11111111111111111111*h0*h1;
1700        const double w60 = -0.19716878364870322056*h1;
1701        const double w61 = -0.19716878364870322056*h0;
1702        const double w62 = -0.052831216351296779436*h0;
1703        const double w63 = -0.052831216351296779436*h1;
1704        const double w64 = 0.19716878364870322056*h1;
1705        const double w65 = 0.052831216351296779436*h1;
1706        const double w66 = 0.19716878364870322056*h0;
1707        const double w67 = 0.052831216351296779436*h0;
1708        const double w68 = -0.5*h1;
1709        const double w69 = -0.5*h0;
1710        const double w70 = 0.5*h1;
1711        const double w71 = 0.5*h0;
1712        const double w72 = 0.1555021169820365539*h0*h1;
1713        const double w73 = 0.041666666666666666667*h0*h1;
1714        const double w74 = 0.01116454968463011277*h0*h1;
1715        const double w75 = 0.25*h0*h1;
1716    
1717        rhs.requireWrite();
1718    #pragma omp parallel
1719        {
1720            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1721    #pragma omp for
1722                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
1723                    for (index_t k0=0; k0<m_NE0; ++k0)  {
1724                        bool add_EM_S=false;
1725                        bool add_EM_F=false;
1726                        vector<double> EM_S(4*4, 0);
1727                        vector<double> EM_F(4, 0);
1728                        const index_t e = k0 + m_NE0*k1;
1729                        ///////////////
1730                        // process A //
1731                        ///////////////
1732                        if (!A.isEmpty()) {
1733                            add_EM_S=true;
1734                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1735                            if (A.actsExpanded()) {
1736                                const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1737                                const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1738                                const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1739                                const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1740                                const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1741                                const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1742                                const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1743                                const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1744                                const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1745                                const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1746                                const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1747                                const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1748                                const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1749                                const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1750                                const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1751                                const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1752                                const double tmp0_0 = A_01_0 + A_01_3;
1753                                const double tmp1_0 = A_00_0 + A_00_1;
1754                                const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1755                                const double tmp3_0 = A_00_2 + A_00_3;
1756                                const double tmp4_0 = A_10_1 + A_10_2;
1757                                const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1758                                const double tmp6_0 = A_01_3 + A_10_0;
1759                                const double tmp7_0 = A_01_0 + A_10_3;
1760                                const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1761                                const double tmp9_0 = A_01_0 + A_10_0;
1762                                const double tmp12_0 = A_11_0 + A_11_2;
1763                                const double tmp10_0 = A_01_3 + A_10_3;
1764                                const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
1765                                const double tmp13_0 = A_01_2 + A_10_1;
1766                                const double tmp11_0 = A_11_1 + A_11_3;
1767                                const double tmp18_0 = A_01_1 + A_10_1;
1768                                const double tmp15_0 = A_01_1 + A_10_2;
1769                                const double tmp16_0 = A_10_0 + A_10_3;
1770                                const double tmp17_0 = A_01_1 + A_01_2;
1771                                const double tmp19_0 = A_01_2 + A_10_2;
1772                                const double tmp0_1 = A_10_3*w8;
1773                                const double tmp1_1 = tmp0_0*w1;
1774                                const double tmp2_1 = A_01_1*w4;
1775                                const double tmp3_1 = tmp1_0*w0;
1776                                const double tmp4_1 = A_01_2*w7;
1777                                const double tmp5_1 = tmp2_0*w3;
1778                                const double tmp6_1 = tmp3_0*w6;
1779                                const double tmp7_1 = A_10_0*w2;
1780                                const double tmp8_1 = tmp4_0*w5;
1781                                const double tmp9_1 = tmp2_0*w10;
1782                                const double tmp14_1 = A_10_0*w8;
1783                                const double tmp23_1 = tmp3_0*w14;
1784                                const double tmp35_1 = A_01_0*w8;
1785                                const double tmp54_1 = tmp13_0*w8;
1786                                const double tmp20_1 = tmp9_0*w4;
1787                                const double tmp25_1 = tmp12_0*w12;
1788                                const double tmp44_1 = tmp7_0*w7;
1789                                const double tmp26_1 = tmp10_0*w4;
1790                                const double tmp52_1 = tmp18_0*w8;
1791                                const double tmp48_1 = A_10_1*w7;
1792                                const double tmp46_1 = A_01_3*w8;
1793                                const double tmp50_1 = A_01_0*w2;
1794                                const double tmp56_1 = tmp19_0*w8;
1795                                const double tmp19_1 = A_10_3*w2;
1796                                const double tmp47_1 = A_10_2*w4;
1797                                const double tmp16_1 = tmp3_0*w0;
1798                                const double tmp18_1 = tmp1_0*w6;
1799                                const double tmp31_1 = tmp11_0*w12;
1800                                const double tmp55_1 = tmp15_0*w2;
1801                                const double tmp39_1 = A_10_2*w7;
1802                                const double tmp11_1 = tmp6_0*w7;
1803                                const double tmp40_1 = tmp11_0*w17;
1804                                const double tmp34_1 = tmp15_0*w8;
1805                                const double tmp33_1 = tmp14_0*w5;
1806                                const double tmp24_1 = tmp11_0*w13;
1807                                const double tmp43_1 = tmp17_0*w5;
1808                                const double tmp15_1 = A_01_2*w4;
1809                                const double tmp53_1 = tmp19_0*w2;
1810                                const double tmp27_1 = tmp3_0*w11;
1811                                const double tmp32_1 = tmp13_0*w2;
1812                                const double tmp10_1 = tmp5_0*w9;
1813                                const double tmp37_1 = A_10_1*w4;
1814                                const double tmp38_1 = tmp5_0*w15;
1815                                const double tmp17_1 = A_01_1*w7;
1816                                const double tmp12_1 = tmp7_0*w4;
1817                                const double tmp22_1 = tmp10_0*w7;
1818                                const double tmp57_1 = tmp18_0*w2;
1819                                const double tmp28_1 = tmp9_0*w7;
1820                                const double tmp29_1 = tmp1_0*w14;
1821                                const double tmp51_1 = tmp11_0*w16;
1822                                const double tmp42_1 = tmp12_0*w16;
1823                                const double tmp49_1 = tmp12_0*w17;
1824                                const double tmp21_1 = tmp1_0*w11;
1825                                const double tmp45_1 = tmp6_0*w4;
1826                                const double tmp13_1 = tmp8_0*w1;
1827                                const double tmp36_1 = tmp16_0*w1;
1828                                const double tmp41_1 = A_01_3*w2;
1829                                const double tmp30_1 = tmp12_0*w13;
1830                                EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1831                                EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
1832                                EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
1833                                EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
1834                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
1835                                EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
1836                                EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
1837                                EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
1838                                EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1839                                EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
1840                                EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
1841                                EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
1842                                EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
1843                                EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
1844                                EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
1845                                EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1846                            } else { // constant data
1847                                const double A_00 = A_p[INDEX2(0,0,2)];
1848                                const double A_10 = A_p[INDEX2(1,0,2)];
1849                                const double A_01 = A_p[INDEX2(0,1,2)];
1850                                const double A_11 = A_p[INDEX2(1,1,2)];
1851                                const double tmp0_0 = A_01 + A_10;
1852                                const double tmp0_1 = A_00*w18;
1853                                const double tmp1_1 = A_01*w19;
1854                                const double tmp2_1 = A_10*w20;
1855                                const double tmp3_1 = A_11*w21;
1856                                const double tmp4_1 = A_00*w22;
1857                                const double tmp5_1 = tmp0_0*w19;
1858                                const double tmp6_1 = A_11*w23;
1859                                const double tmp7_1 = A_11*w25;
1860                                const double tmp8_1 = A_00*w24;
1861                                const double tmp9_1 = tmp0_0*w20;
1862                                const double tmp10_1 = A_01*w20;
1863                                const double tmp11_1 = A_11*w27;
1864                                const double tmp12_1 = A_00*w26;
1865                                const double tmp13_1 = A_10*w19;
1866                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1867                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1868                                EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1869                                EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1870                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1871                                EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1872                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1873                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1874                                EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1875                                EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1876                                EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1877                                EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1878                                EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1879                                EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1880                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1881                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1882                            }
1883                        }
1884                        ///////////////
1885                        // process B //
1886                        ///////////////
1887                        if (!B.isEmpty()) {
1888                            add_EM_S=true;
1889                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1890                            if (B.actsExpanded()) {
1891                                const double B_0_0 = B_p[INDEX2(0,0,2)];
1892                                const double B_1_0 = B_p[INDEX2(1,0,2)];
1893                                const double B_0_1 = B_p[INDEX2(0,1,2)];
1894                                const double B_1_1 = B_p[INDEX2(1,1,2)];
1895                                const double B_0_2 = B_p[INDEX2(0,2,2)];
1896                                const double B_1_2 = B_p[INDEX2(1,2,2)];
1897                                const double B_0_3 = B_p[INDEX2(0,3,2)];
1898                                const double B_1_3 = B_p[INDEX2(1,3,2)];
1899                                const double tmp0_0 = B_1_0 + B_1_1;
1900                                const double tmp1_0 = B_1_2 + B_1_3;
1901                                const double tmp2_0 = B_0_1 + B_0_3;
1902                                const double tmp3_0 = B_0_0 + B_0_2;
1903                                const double tmp63_1 = B_1_1*w42;
1904                                const double tmp79_1 = B_1_1*w40;
1905                                const double tmp37_1 = tmp3_0*w35;
1906                                const double tmp8_1 = tmp0_0*w32;
1907                                const double tmp71_1 = B_0_1*w34;
1908                                const double tmp19_1 = B_0_3*w31;
1909                                const double tmp15_1 = B_0_3*w34;
1910                                const double tmp9_1 = tmp3_0*w34;
1911                                const double tmp35_1 = B_1_0*w36;
1912                                const double tmp66_1 = B_0_3*w28;
1913                                const double tmp28_1 = B_1_0*w42;
1914                                const double tmp22_1 = B_1_0*w40;
1915                                const double tmp16_1 = B_1_2*w29;
1916                                const double tmp6_1 = tmp2_0*w35;
1917                                const double tmp55_1 = B_1_3*w40;
1918                                const double tmp50_1 = B_1_3*w42;
1919                                const double tmp7_1 = tmp1_0*w29;
1920                                const double tmp1_1 = tmp1_0*w32;
1921                                const double tmp57_1 = B_0_3*w30;
1922                                const double tmp18_1 = B_1_1*w32;
1923                                const double tmp53_1 = B_1_0*w41;
1924                                const double tmp61_1 = B_1_3*w36;
1925                                const double tmp27_1 = B_0_3*w38;
1926                                const double tmp64_1 = B_0_2*w30;
1927                                const double tmp76_1 = B_0_1*w38;
1928                                const double tmp39_1 = tmp2_0*w34;
1929                                const double tmp62_1 = B_0_1*w31;
1930                                const double tmp56_1 = B_0_0*w31;
1931                                const double tmp49_1 = B_1_1*w36;
1932                                const double tmp2_1 = B_0_2*w31;
1933                                const double tmp23_1 = B_0_2*w33;
1934                                const double tmp38_1 = B_1_1*w43;
1935                                const double tmp74_1 = B_1_2*w41;
1936                                const double tmp43_1 = B_1_1*w41;
1937                                const double tmp58_1 = B_0_2*w28;
1938                                const double tmp67_1 = B_0_0*w33;
1939                                const double tmp33_1 = tmp0_0*w39;
1940                                const double tmp4_1 = B_0_0*w28;
1941                                const double tmp20_1 = B_0_0*w30;
1942                                const double tmp13_1 = B_0_2*w38;
1943                                const double tmp65_1 = B_1_2*w43;
1944                                const double tmp0_1 = tmp0_0*w29;
1945                                const double tmp41_1 = tmp3_0*w33;
1946                                const double tmp73_1 = B_0_2*w37;
1947                                const double tmp69_1 = B_0_0*w38;
1948                                const double tmp48_1 = B_1_2*w39;
1949                                const double tmp59_1 = B_0_1*w33;
1950                                const double tmp17_1 = B_1_3*w41;
1951                                const double tmp5_1 = B_0_3*w33;
1952                                const double tmp3_1 = B_0_1*w30;
1953                                const double tmp21_1 = B_0_1*w28;
1954                                const double tmp42_1 = B_1_0*w29;
1955                                const double tmp54_1 = B_1_2*w32;
1956                                const double tmp60_1 = B_1_0*w39;
1957                                const double tmp32_1 = tmp1_0*w36;
1958                                const double tmp10_1 = B_0_1*w37;
1959                                const double tmp14_1 = B_0_0*w35;
1960                                const double tmp29_1 = B_0_1*w35;
1961                                const double tmp26_1 = B_1_2*w36;
1962                                const double tmp30_1 = B_1_3*w43;
1963                                const double tmp70_1 = B_0_2*w35;
1964                                const double tmp34_1 = B_1_3*w39;
1965                                const double tmp51_1 = B_1_0*w43;
1966                                const double tmp31_1 = B_0_2*w34;
1967                                const double tmp45_1 = tmp3_0*w28;
1968                                const double tmp11_1 = tmp1_0*w39;
1969                                const double tmp52_1 = B_1_1*w29;
1970                                const double tmp44_1 = B_1_3*w32;
1971                                const double tmp25_1 = B_1_1*w39;
1972                                const double tmp47_1 = tmp2_0*w33;
1973                                const double tmp72_1 = B_1_3*w29;
1974                                const double tmp40_1 = tmp2_0*w28;
1975                                const double tmp46_1 = B_1_2*w40;
1976                                const double tmp36_1 = B_1_2*w42;
1977                                const double tmp24_1 = B_0_0*w37;
1978                                const double tmp77_1 = B_0_3*w35;
1979                                const double tmp68_1 = B_0_3*w37;
1980                                const double tmp78_1 = B_0_0*w34;
1981                                const double tmp12_1 = tmp0_0*w36;
1982                                const double tmp75_1 = B_1_0*w32;
1983                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1984                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1985                                EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1986                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
1987                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1988                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1989                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1990                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1991                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1992                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1993                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1994                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1995                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
1996                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1997                                EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1998                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1999                            } else { // constant data
2000                                const double B_0 = B_p[0];
2001                                const double B_1 = B_p[1];
2002                                const double tmp0_1 = B_0*w44;
2003                                const double tmp1_1 = B_1*w45;
2004                                const double tmp2_1 = B_0*w46;
2005                                const double tmp3_1 = B_0*w47;
2006                                const double tmp4_1 = B_1*w48;
2007                                const double tmp5_1 = B_1*w49;
2008                                const double tmp6_1 = B_1*w50;
2009                                const double tmp7_1 = B_0*w51;
2010                                EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
2011                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2012                                EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
2013                                EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
2014                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2015                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
2016                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
2017                                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
2018                                EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
2019                                EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
2020                                EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
2021                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
2022                                EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
2023                                EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
2024                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
2025                                EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
2026                            }
2027                        }
2028                        ///////////////
2029                        // process C //
2030                        ///////////////
2031                        if (!C.isEmpty()) {
2032                            add_EM_S=true;
2033                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2034                            if (C.actsExpanded()) {
2035                                const double C_0_0 = C_p[INDEX2(0,0,2)];
2036                                const double C_1_0 = C_p[INDEX2(1,0,2)];
2037                                const double C_0_1 = C_p[INDEX2(0,1,2)];
2038                                const double C_1_1 = C_p[INDEX2(1,1,2)];
2039                                const double C_0_2 = C_p[INDEX2(0,2,2)];
2040                                const double C_1_2 = C_p[INDEX2(1,2,2)];
2041                                const double C_0_3 = C_p[INDEX2(0,3,2)];
2042                                const double C_1_3 = C_p[INDEX2(1,3,2)];
2043                                const double tmp0_0 = C_1_0 + C_1_1;
2044                                const double tmp1_0 = C_1_2 + C_1_3;
2045                                const double tmp2_0 = C_0_1 + C_0_3;
2046                                const double tmp3_0 = C_0_0 + C_0_2;
2047                                const double tmp64_1 = C_0_2*w30;
2048                                const double tmp14_1 = C_0_2*w28;
2049                                const double tmp19_1 = C_0_3*w31;
2050                                const double tmp22_1 = C_1_0*w40;
2051                                const double tmp37_1 = tmp3_0*w35;
2052                                const double tmp29_1 = C_0_1*w35;
2053                                const double tmp73_1 = C_0_2*w37;
2054                                const double tmp74_1 = C_1_2*w41;
2055                                const double tmp52_1 = C_1_3*w39;
2056                                const double tmp25_1 = C_1_1*w39;
2057                                const double tmp62_1 = C_0_1*w31;
2058                                const double tmp79_1 = C_1_1*w40;
2059                                const double tmp43_1 = C_1_1*w36;
2060                                const double tmp27_1 = C_0_3*w38;
2061                                const double tmp28_1 = C_1_0*w42;
2062                                const double tmp63_1 = C_1_1*w42;
2063                                const double tmp59_1 = C_0_3*w34;
2064                                const double tmp72_1 = C_1_3*w29;
2065                                const double tmp40_1 = tmp2_0*w35;
2066                                const double tmp13_1 = C_0_3*w30;
2067                                const double tmp51_1 = C_1_2*w40;
2068                                const double tmp54_1 = C_1_2*w42;
2069                                const double tmp12_1 = C_0_0*w31;
2070                                const double tmp2_1 = tmp1_0*w32;
2071                                const double tmp68_1 = C_0_2*w31;
2072                                const double tmp75_1 = C_1_0*w32;
2073                                const double tmp49_1 = C_1_1*w41;
2074                                const double tmp4_1 = C_0_2*w35;
2075                                const double tmp66_1 = C_0_3*w28;
2076                                const double tmp56_1 = C_0_1*w37;
2077                                const double tmp5_1 = C_0_1*w34;
2078                                const double tmp38_1 = tmp2_0*w34;
2079                                const double tmp76_1 = C_0_1*w38;
2080                                const double tmp21_1 = C_0_1*w28;
2081                                const double tmp69_1 = C_0_1*w30;
2082                                const double tmp53_1 = C_1_0*w36;
2083                                const double tmp42_1 = C_1_2*w39;
2084                                const double tmp32_1 = tmp1_0*w29;
2085                                const double tmp45_1 = C_1_0*w43;
2086                                const double tmp33_1 = tmp0_0*w32;
2087                                const double tmp35_1 = C_1_0*w41;
2088                                const double tmp26_1 = C_1_2*w36;
2089                                const double tmp67_1 = C_0_0*w33;
2090                                const double tmp31_1 = C_0_2*w34;
2091                                const double tmp20_1 = C_0_0*w30;
2092                                const double tmp70_1 = C_0_0*w28;
2093                                const double tmp8_1 = tmp0_0*w39;
2094                                const double tmp30_1 = C_1_3*w43;
2095                                const double tmp0_1 = tmp0_0*w29;
2096                                const double tmp17_1 = C_1_3*w41;
2097                                const double tmp58_1 = C_0_0*w35;
2098                                const double tmp9_1 = tmp3_0*w33;
2099                                const double tmp61_1 = C_1_3*w36;
2100                                const double tmp41_1 = tmp3_0*w34;
2101                                const double tmp50_1 = C_1_3*w32;
2102                                const double tmp18_1 = C_1_1*w32;
2103                                const double tmp6_1 = tmp1_0*w36;
2104                                const double tmp3_1 = C_0_0*w38;
2105                                const double tmp34_1 = C_1_1*w29;
2106                                const double tmp77_1 = C_0_3*w35;
2107                                const double tmp65_1 = C_1_2*w43;
2108                                const double tmp71_1 = C_0_3*w33;
2109                                const double tmp55_1 = C_1_1*w43;
2110                                const double tmp46_1 = tmp3_0*w28;
2111                                const double tmp24_1 = C_0_0*w37;
2112                                const double tmp10_1 = tmp1_0*w39;
2113                                const double tmp48_1 = C_1_0*w29;
2114                                const double tmp15_1 = C_0_1*w33;
2115                                const double tmp36_1 = C_1_2*w32;
2116                                const double tmp60_1 = C_1_0*w39;
2117                                const double tmp47_1 = tmp2_0*w33;
2118                                const double tmp16_1 = C_1_2*w29;
2119                                const double tmp1_1 = C_0_3*w37;
2120                                const double tmp7_1 = tmp2_0*w28;
2121                                const double tmp39_1 = C_1_3*w40;
2122                                const double tmp44_1 = C_1_3*w42;
2123                                const double tmp57_1 = C_0_2*w38;
2124                                const double tmp78_1 = C_0_0*w34;
2125                                const double tmp11_1 = tmp0_0*w36;
2126                                const double tmp23_1 = C_0_2*w33;
2127                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
2128                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2129                                EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2130                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
2131                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2132                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
2133                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
2134                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
2135                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
2136                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2137                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
2138                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2139                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
2140                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
2141                                EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2142                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2143                            } else { // constant data
2144                                const double C_0 = C_p[0];
2145                                const double C_1 = C_p[1];
2146                                const double tmp0_1 = C_0*w47;
2147                                const double tmp1_1 = C_1*w45;
2148                                const double tmp2_1 = C_1*w48;
2149                                const double tmp3_1 = C_0*w51;
2150                                const double tmp4_1 = C_0*w44;
2151                                const double tmp5_1 = C_1*w49;
2152                                const double tmp6_1 = C_1*w50;
2153                                const double tmp7_1 = C_0*w46;
2154                                EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
2155                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
2156                                EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
2157                                EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
2158                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2159                                EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
2160                                EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
2161                                EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
2162                                EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
2163                                EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
2164                                EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
2165                                EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
2166                                EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
2167                                EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
2168                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2169                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
2170                            }
2171                        }
2172                        ///////////////
2173                        // process D //
2174                        ///////////////
2175                        if (!D.isEmpty()) {
2176                            add_EM_S=true;
2177                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2178                            if (D.actsExpanded()) {
2179                                const double D_0 = D_p[0];
2180                                const double D_1 = D_p[1];
2181                                const double D_2 = D_p[2];
2182                                const double D_3 = D_p[3];
2183                                const double tmp0_0 = D_0 + D_1;
2184                                const double tmp1_0 = D_2 + D_3;
2185                                const double tmp2_0 = D_0 + D_1 + D_2 + D_3;
2186                                const double tmp3_0 = D_1 + D_2;
2187                                const double tmp4_0 = D_1 + D_3;
2188                                const double tmp5_0 = D_0 + D_2;
2189                                const double tmp6_0 = D_0 + D_3;
2190                                const double tmp0_1 = tmp0_0*w52;
2191                                const double tmp1_1 = tmp1_0*w53;
2192                                const double tmp2_1 = tmp2_0*w54;
2193                                const double tmp3_1 = tmp1_0*w52;
2194                                const double tmp4_1 = tmp0_0*w53;
2195                                const double tmp5_1 = tmp3_0*w54;
2196                                const double tmp6_1 = D_0*w55;
2197                                const double tmp7_1 = D_3*w56;
2198                                const double tmp8_1 = D_3*w55;
2199                                const double tmp9_1 = D_0*w56;
2200                                const double tmp10_1 = tmp4_0*w52;
2201                                const double tmp11_1 = tmp5_0*w53;
2202                                const double tmp12_1 = tmp5_0*w52;
2203                                const double tmp13_1 = tmp4_0*w53;
2204                                const double tmp14_1 = tmp6_0*w54;
2205                                const double tmp15_1 = D_2*w55;
2206                                const double tmp16_1 = D_1*w56;
2207                                const double tmp17_1 = D_1*w55;
2208                                const double tmp18_1 = D_2*w56;
2209                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2210                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2211                                EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
2212                                EM_S[INDEX2(3,0,4)]+=tmp2_1;
2213                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2214                                EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
2215                                EM_S[INDEX2(2,1,4)]+=tmp2_1;
2216                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
2217                                EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
2218                                EM_S[INDEX2(1,2,4)]+=tmp2_1;
2219                                EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
2220                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
2221                                EM_S[INDEX2(0,3,4)]+=tmp2_1;
2222                                EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
2223                                EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2224                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2225                            } else { // constant data
2226                                const double tmp0_1 = D_p[0]*w57;
2227                                const double tmp1_1 = D_p[0]*w58;
2228                                const double tmp2_1 = D_p[0]*w59;
2229                                EM_S[INDEX2(0,0,4)]+=tmp2_1;
2230                                EM_S[INDEX2(1,0,4)]+=tmp0_1;
2231                                EM_S[INDEX2(2,0,4)]+=tmp0_1;
2232                                EM_S[INDEX2(3,0,4)]+=tmp1_1;
2233                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
2234                                EM_S[INDEX2(1,1,4)]+=tmp2_1;
2235                                EM_S[INDEX2(2,1,4)]+=tmp1_1;
2236                                EM_S[INDEX2(3,1,4)]+=tmp0_1;
2237                                EM_S[INDEX2(0,2,4)]+=tmp0_1;
2238                                EM_S[INDEX2(1,2,4)]+=tmp1_1;
2239                                EM_S[INDEX2(2,2,4)]+=tmp2_1;
2240                                EM_S[INDEX2(3,2,4)]+=tmp0_1;
2241                                EM_S[INDEX2(0,3,4)]+=tmp1_1;
2242                                EM_S[INDEX2(1,3,4)]+=tmp0_1;
2243                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
2244                                EM_S[INDEX2(3,3,4)]+=tmp2_1;
2245                            }
2246                        }
2247                        ///////////////
2248                        // process X //
2249                        ///////////////
2250                        if (!X.isEmpty()) {
2251                            add_EM_F=true;
2252                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2253                            if (X.actsExpanded()) {
2254                                const double X_0_0 = X_p[INDEX2(0,0,2)];
2255                                const double X_1_0 = X_p[INDEX2(1,0,2)];
2256                                const double X_0_1 = X_p[INDEX2(0,1,2)];
2257                                const double X_1_1 = X_p[INDEX2(1,1,2)];
2258                                const double X_0_2 = X_p[INDEX2(0,2,2)];
2259                                const double X_1_2 = X_p[INDEX2(1,2,2)];
2260                                const double X_0_3 = X_p[INDEX2(0,3,2)];
2261                                const double X_1_3 = X_p[INDEX2(1,3,2)];
2262                                const double tmp0_0 = X_0_2 + X_0_3;
2263                                const double tmp1_0 = X_1_1 + X_1_3;
2264                                const double tmp2_0 = X_1_0 + X_1_2;
2265                                const double tmp3_0 = X_0_0 + X_0_1;
2266                                const double tmp0_1 = tmp0_0*w63;
2267                                const double tmp1_1 = tmp1_0*w62;
2268                                const double tmp2_1 = tmp2_0*w61;
2269                                const double tmp3_1 = tmp3_0*w60;
2270                                const double tmp4_1 = tmp0_0*w65;
2271                                const double tmp5_1 = tmp3_0*w64;
2272                                const double tmp6_1 = tmp2_0*w62;
2273                                const double tmp7_1 = tmp1_0*w61;
2274                                const double tmp8_1 = tmp2_0*w66;
2275                                const double tmp9_1 = tmp3_0*w63;
2276                                const double tmp10_1 = tmp0_0*w60;
2277                                const double tmp11_1 = tmp1_0*w67;
2278                                const double tmp12_1 = tmp1_0*w66;
2279                                const double tmp13_1 = tmp3_0*w65;
2280                                const double tmp14_1 = tmp0_0*w64;
2281                                const double tmp15_1 = tmp2_0*w67;
2282                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2283                                EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2284                                EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2285                                EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2286                            } else { // constant data
2287                                const double tmp0_1 = X_p[1]*w69;
2288                                const double tmp1_1 = X_p[0]*w68;
2289                                const double tmp2_1 = X_p[0]*w70;
2290                                const double tmp3_1 = X_p[1]*w71;
2291                                EM_F[0]+=tmp0_1 + tmp1_1;
2292                                EM_F[1]+=tmp0_1 + tmp2_1;
2293                                EM_F[2]+=tmp1_1 + tmp3_1;
2294                                EM_F[3]+=tmp2_1 + tmp3_1;
2295                            }
2296                        }
2297                        ///////////////
2298                        // process Y //
2299                        ///////////////
2300                        if (!Y.isEmpty()) {
2301                            add_EM_F=true;
2302                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2303                            if (Y.actsExpanded()) {
2304                                const double Y_0 = Y_p[0];
2305                                const double Y_1 = Y_p[1];
2306                                const double Y_2 = Y_p[2];
2307                                const double Y_3 = Y_p[3];
2308                                const double tmp0_0 = Y_1 + Y_2;
2309                                const double tmp1_0 = Y_0 + Y_3;
2310                                const double tmp0_1 = Y_0*w72;
2311                                const double tmp1_1 = tmp0_0*w73;
2312                                const double tmp2_1 = Y_3*w74;
2313                                const double tmp3_1 = Y_1*w72;
2314                                const double tmp4_1 = tmp1_0*w73;
2315                                const double tmp5_1 = Y_2*w74;
2316                                const double tmp6_1 = Y_2*w72;
2317                                const double tmp7_1 = Y_1*w74;
2318                                const double tmp8_1 = Y_3*w72;
2319                                const double tmp9_1 = Y_0*w74;
2320                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2321                                EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2322                                EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2323                                EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2324                            } else { // constant data
2325                                const double tmp0_1 = Y_p[0]*w75;
2326                                EM_F[0]+=tmp0_1;
2327                                EM_F[1]+=tmp0_1;
2328                                EM_F[2]+=tmp0_1;
2329                                EM_F[3]+=tmp0_1;
2330                            }
2331                        }
2332    
2333                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2334                        const index_t firstNode=m_N0*k1+k0;
2335                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2336                    } // end k0 loop
2337                } // end k1 loop
2338            } // end of colouring
2339        } // end of parallel region
2340    }
2341    
2342    //protected
2343    void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2344            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2345            const escript::Data& C, const escript::Data& D,
2346            const escript::Data& X, const escript::Data& Y) const
2347    {
2348        const double h0 = m_l0/m_gNE0;
2349        const double h1 = m_l1/m_gNE1;
2350        const double w0 = -.25*h1/h0;
2351        const double w1 = .25;
2352        const double w2 = -.25;
2353        const double w3 = .25*h0/h1;
2354        const double w4 = -.25*h0/h1;
2355        const double w5 = .25*h1/h0;
2356        const double w6 = -.125*h1;
2357        const double w7 = -.125*h0;
2358        const double w8 = .125*h1;
2359        const double w9 = .125*h0;
2360        const double w10 = .0625*h0*h1;
2361        const double w11 = -.5*h1;
2362        const double w12 = -.5*h0;
2363        const double w13 = .5*h1;
2364        const double w14 = .5*h0;
2365        const double w15 = .25*h0*h1;
2366    
2367        rhs.requireWrite();
2368    #pragma omp parallel
2369        {
2370            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2371    #pragma omp for
2372                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2373                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2374                        bool add_EM_S=false;
2375                        bool add_EM_F=false;
2376                        vector<double> EM_S(4*4, 0);
2377                        vector<double> EM_F(4, 0);
2378                        const index_t e = k0 + m_NE0*k1;
2379                        ///////////////
2380                        // process A //
2381                        ///////////////
2382                        if (!A.isEmpty()) {
2383                            add_EM_S=true;
2384                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2385                            const double A_00 = A_p[INDEX2(0,0,2)];
2386                            const double A_10 = A_p[INDEX2(1,0,2)];
2387                            const double A_01 = A_p[INDEX2(0,1,2)];
2388                            const double A_11 = A_p[INDEX2(1,1,2)];
2389                            const double tmp0_0 = A_01 + A_10;
2390                            const double tmp0_1 = A_11*w3;
2391                            const double tmp1_1 = A_00*w0;
2392                            const double tmp2_1 = A_01*w1;
2393                            const double tmp3_1 = A_10*w2;
2394                            const double tmp4_1 = tmp0_0*w1;
2395                            const double tmp5_1 = A_11*w4;
2396                            const double tmp6_1 = A_00*w5;
2397                            const double tmp7_1 = tmp0_0*w2;
2398                            const double tmp8_1 = A_10*w1;
2399                            const double tmp9_1 = A_01*w2;
2400                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2401                            EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2402                            EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2403                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2404                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2405                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2406                            EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2407                            EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2408                            EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2409                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2410                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2411                            EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2412                            EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2413                            EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2414                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2415                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2416                        }
2417                        ///////////////
2418                        // process B //
2419                        ///////////////
2420                        if (!B.isEmpty()) {
2421                            add_EM_S=true;
2422                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2423                            const double tmp2_1 = B_p[0]*w8;
2424                            const double tmp0_1 = B_p[0]*w6;
2425                            const double tmp3_1 = B_p[1]*w9;
2426                            const double tmp1_1 = B_p[1]*w7;
2427                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2428                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2429                            EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
2430                            EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;
2431                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2432                            EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;
2433                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;
2434                            EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;
2435                            EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;
2436                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
2437                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;
2438                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2439                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;
2440                            EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;
2441                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;
2442                            EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
2443                        }
2444                        ///////////////
2445                        // process C //
2446                        ///////////////
2447                        if (!C.isEmpty()) {
2448                            add_EM_S=true;
2449                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2450                            const double tmp1_1 = C_p[1]*w7;
2451                            const double tmp0_1 = C_p[0]*w8;
2452                            const double tmp3_1 = C_p[0]*w6;
2453                            const double tmp2_1 = C_p[1]*w9;
2454                            EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2455                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2456                            EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
2457                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
2458                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2459                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;
2460                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;
2461                            EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;
2462                            EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;
2463                            EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
2464                            EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
2465                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2466                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;
2467                            EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;
2468                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2469                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;
2470                        }
2471                        ///////////////
2472                        // process D //
2473                        ///////////////
2474                        if (!D.isEmpty()) {
2475                            add_EM_S=true;
2476                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2477                            const double tmp0_1 = D_p[0]*w10;
2478                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
2479                            EM_S[INDEX2(1,0,4)]+=tmp0_1;
2480                            EM_S[INDEX2(2,0,4)]+=tmp0_1;
2481                            EM_S[INDEX2(3,0,4)]+=tmp0_1;
2482                            EM_S[INDEX2(0,1,4)]+=tmp0_1;
2483                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
2484                            EM_S[INDEX2(2,1,4)]+=tmp0_1;
2485                            EM_S[INDEX2(3,1,4)]+=tmp0_1;
2486                            EM_S[INDEX2(0,2,4)]+=tmp0_1;
2487                            EM_S[INDEX2(1,2,4)]+=tmp0_1;
2488                            EM_S[INDEX2(2,2,4)]+=tmp0_1;
2489                            EM_S[INDEX2(3,2,4)]+=tmp0_1;
2490                            EM_S[INDEX2(0,3,4)]+=tmp0_1;
2491                            EM_S[INDEX2(1,3,4)]+=tmp0_1;
2492                            EM_S[INDEX2(2,3,4)]+=tmp0_1;
2493                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
2494                        }
2495                        ///////////////
2496                        // process X //
2497                        ///////////////
2498                        if (!X.isEmpty()) {
2499                            add_EM_F=true;
2500                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2501                            const double tmp0_1 = X_p[0]*w11;
2502                            const double tmp2_1 = X_p[0]*w13;
2503                            const double tmp1_1 = X_p[1]*w12;
2504                            const double tmp3_1 = X_p[1]*w14;
2505                            EM_F[0]+=tmp0_1 + tmp1_1;
2506                            EM_F[1]+=tmp1_1 + tmp2_1;
2507                            EM_F[2]+=tmp0_1 + tmp3_1;
2508                            EM_F[3]+=tmp2_1 + tmp3_1;
2509                        }
2510                        ///////////////
2511                        // process Y //
2512                        ///////////////
2513                        if (!Y.isEmpty()) {
2514                            add_EM_F=true;
2515                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2516                            const double tmp0_1 = Y_p[0]*w15;
2517                            EM_F[0]+=tmp0_1;
2518                            EM_F[1]+=tmp0_1;
2519                            EM_F[2]+=tmp0_1;
2520                            EM_F[3]+=tmp0_1;
2521                        }
2522    
2523                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2524                        const index_t firstNode=m_N0*k1+k0;
2525                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2526                    } // end k0 loop
2527                } // end k1 loop
2528            } // end of colouring
2529        } // end of parallel region
2530    }
2531    
2532    //protected
2533    void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2534            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2535            const escript::Data& C, const escript::Data& D,
2536            const escript::Data& X, const escript::Data& Y) const
2537    {
2538        const double h0 = m_l0/m_gNE0;
2539        const double h1 = m_l1/m_gNE1;
2540        dim_t numEq, numComp;
2541        if (!mat)
2542            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
2543        else {
2544            numEq=mat->logical_row_block_size;
2545            numComp=mat->logical_col_block_size;
2546        }
2547    
2548        const double w0 = -0.1555021169820365539*h1/h0;
2549        const double w1 = 0.041666666666666666667;
2550        const double w2 = -0.15550211698203655390;
2551        const double w3 = 0.041666666666666666667*h0/h1;
2552        const double w4 = 0.15550211698203655390;
2553        const double w5 = -0.041666666666666666667;
2554        const double w6 = -0.01116454968463011277*h1/h0;
2555        const double w7 = 0.011164549684630112770;
2556        const double w8 = -0.011164549684630112770;
2557        const double w9 = -0.041666666666666666667*h1/h0;
2558        const double w10 = -0.041666666666666666667*h0/h1;
2559        const double w11 = 0.1555021169820365539*h1/h0;
2560        const double w12 = 0.1555021169820365539*h0/h1;
2561        const double w13 = 0.01116454968463011277*h0/h1;
2562        const double w14 = 0.01116454968463011277*h1/h0;
2563        const double w15 = 0.041666666666666666667*h1/h0;
2564        const double w16 = -0.01116454968463011277*h0/h1;
2565        const double w17 = -0.1555021169820365539*h0/h1;
2566        const double w18 = -0.33333333333333333333*h1/h0;
2567        const double w19 = 0.25000000000000000000;
2568        const double w20 = -0.25000000000000000000;
2569        const double w21 = 0.16666666666666666667*h0/h1;
2570        const double w22 = -0.16666666666666666667*h1/h0;
2571        const double w23 = -0.16666666666666666667*h0/h1;
2572        const double w24 = 0.33333333333333333333*h1/h0;
2573        const double w25 = 0.33333333333333333333*h0/h1;
2574        const double w26 = 0.16666666666666666667*h1/h0;
2575        const double w27 = -0.33333333333333333333*h0/h1;
2576        const double w28 = -0.032861463941450536761*h1;
2577        const double w29 = -0.032861463941450536761*h0;
2578        const double w30 = -0.12264065304058601714*h1;
2579        const double w31 = -0.0023593469594139828636*h1;
2580        const double w32 = -0.008805202725216129906*h0;
2581        const double w33 = -0.008805202725216129906*h1;
2582        const double w34 = 0.032861463941450536761*h1;
2583        const double w35 = 0.008805202725216129906*h1;
2584        const double w36 = 0.008805202725216129906*h0;
2585        const double w37 = 0.0023593469594139828636*h1;
2586        const double w38 = 0.12264065304058601714*h1;
2587        const double w39 = 0.032861463941450536761*h0;
2588        const double w40 = -0.12264065304058601714*h0;
2589        const double w41 = -0.0023593469594139828636*h0;
2590        const double w42 = 0.0023593469594139828636*h0;
2591        const double w43 = 0.12264065304058601714*h0;
2592        const double w44 = -0.16666666666666666667*h1;
2593        const double w45 = -0.083333333333333333333*h0;
2594        const double w46 = 0.083333333333333333333*h1;
2595        const double w47 = 0.16666666666666666667*h1;
2596        const double w48 = 0.083333333333333333333*h0;
2597        const double w49 = -0.16666666666666666667*h0;
2598        const double w50 = 0.16666666666666666667*h0;
2599        const double w51 = -0.083333333333333333333*h1;
2600        const double w52 = 0.025917019497006092316*h0*h1;
2601        const double w53 = 0.0018607582807716854616*h0*h1;
2602        const double w54 = 0.0069444444444444444444*h0*h1;
2603        const double w55 = 0.09672363354357992482*h0*h1;
2604        const double w56 = 0.00049858867864229740201*h0*h1;
2605        const double w57 = 0.055555555555555555556*h0*h1;
2606        const double w58 = 0.027777777777777777778*h0*h1;
2607        const double w59 = 0.11111111111111111111*h0*h1;
2608        const double w60 = -0.19716878364870322056*h1;
2609        const double w61 = -0.19716878364870322056*h0;
2610        const double w62 = -0.052831216351296779436*h0;
2611        const double w63 = -0.052831216351296779436*h1;
2612        const double w64 = 0.19716878364870322056*h1;
2613        const double w65 = 0.052831216351296779436*h1;
2614        const double w66 = 0.19716878364870322056*h0;
2615        const double w67 = 0.052831216351296779436*h0;
2616        const double w68 = -0.5*h1;
2617        const double w69 = -0.5*h0;
2618        const double w70 = 0.5*h1;
2619        const double w71 = 0.5*h0;
2620        const double w72 = 0.1555021169820365539*h0*h1;
2621        const double w73 = 0.041666666666666666667*h0*h1;
2622        const double w74 = 0.01116454968463011277*h0*h1;
2623        const double w75 = 0.25*h0*h1;
2624    
2625        rhs.requireWrite();
2626    #pragma omp parallel
2627        {
2628            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2629    #pragma omp for
2630                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2631                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2632                        bool add_EM_S=false;
2633                        bool add_EM_F=false;
2634                        vector<double> EM_S(4*4*numEq*numComp, 0);
2635                        vector<double> EM_F(4*numEq, 0);
2636                        const index_t e = k0 + m_NE0*k1;
2637                        ///////////////
2638                        // process A //
2639                        ///////////////
2640                        if (!A.isEmpty()) {
2641                            add_EM_S=true;
2642                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2643                            if (A.actsExpanded()) {
2644                                for (index_t k=0; k<numEq; k++) {
2645                                    for (index_t m=0; m<numComp; m++) {
2646                                        const double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];
2647                                        const double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];
2648                                        const double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];
2649                                        const double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];
2650                                        const double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];
2651                                        const double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];
2652                                        const double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];
2653                                        const double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];
2654                                        const double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];
2655                                        const double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];
2656                                        const double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];
2657                                        const double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];
2658                                        const double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];
2659                                        const double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];
2660                                        const double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];
2661                                        const double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];
2662                                        const double tmp0_0 = A_01_0 + A_01_3;
2663                                        const double tmp1_0 = A_00_0 + A_00_1;
2664                                        const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
2665                                        const double tmp3_0 = A_00_2 + A_00_3;
2666                                        const double tmp4_0 = A_10_1 + A_10_2;
2667                                        const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
2668                                        const double tmp6_0 = A_01_3 + A_10_0;
2669                                        const double tmp7_0 = A_01_0 + A_10_3;
2670                                        const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
2671                                        const double tmp9_0 = A_01_0 + A_10_0;
2672                                        const double tmp10_0 = A_01_3 + A_10_3;
2673                                        const double tmp11_0 = A_11_1 + A_11_3;
2674                                        const double tmp12_0 = A_11_0 + A_11_2;
2675                                        const double tmp13_0 = A_01_2 + A_10_1;
2676                                        const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
2677                                        const double tmp15_0 = A_01_1 + A_10_2;
2678                                        const double tmp16_0 = A_10_0 + A_10_3;
2679                                        const double tmp17_0 = A_01_1 + A_01_2;
2680                                        const double tmp18_0 = A_01_1 + A_10_1;
2681                                        const double tmp19_0 = A_01_2 + A_10_2;
2682                                        const double tmp0_1 = A_10_3*w8;
2683                                        const double tmp1_1 = tmp0_0*w1;
2684                                        const double tmp2_1 = A_01_1*w4;
2685                                        const double tmp3_1 = tmp1_0*w0;
2686                                        const double tmp4_1 = A_01_2*w7;
2687                                        const double tmp5_1 = tmp2_0*w3;
2688                                        const double tmp6_1 = tmp3_0*w6;
2689                                        const double tmp7_1 = A_10_0*w2;
2690                                        const double tmp8_1 = tmp4_0*w5;
2691                                        const double tmp9_1 = tmp2_0*w10;
2692                                        const double tmp10_1 = tmp5_0*w9;
2693                                        const double tmp11_1 = tmp6_0*w7;
2694                                        const double tmp12_1 = tmp7_0*w4;
2695                                        const double tmp13_1 = tmp8_0*w1;
2696                                        const double tmp14_1 = A_10_0*w8;
2697                                        const double tmp15_1 = A_01_2*w4;
2698                                        const double tmp16_1 = tmp3_0*w0;
2699                                        const double tmp17_1 = A_01_1*w7;
2700                                        const double tmp18_1 = tmp1_0*w6;
2701                                        const double tmp19_1 = A_10_3*w2;
2702                                        const double tmp20_1 = tmp9_0*w4;
2703                                        const double tmp21_1 = tmp1_0*w11;
2704                                        const double tmp22_1 = tmp10_0*w7;
2705                                        const double tmp23_1 = tmp3_0*w14;
2706                                        const double tmp24_1 = tmp11_0*w13;
2707                                        const double tmp25_1 = tmp12_0*w12;
2708                                        const double tmp26_1 = tmp10_0*w4;
2709                                        const double tmp27_1 = tmp3_0*w11;
2710                                        const double tmp28_1 = tmp9_0*w7;
2711                                        const double tmp29_1 = tmp1_0*w14;
2712                                        const double tmp30_1 = tmp12_0*w13;
2713                                        const double tmp31_1 = tmp11_0*w12;
2714                                        const double tmp32_1 = tmp13_0*w2;
2715                                        const double tmp33_1 = tmp14_0*w5;
2716                                        const double tmp34_1 = tmp15_0*w8;
2717                                        const double tmp35_1 = A_01_0*w8;
2718                                        const double tmp36_1 = tmp16_0*w1;
2719                                        const double tmp37_1 = A_10_1*w4;
2720                                        const double tmp38_1 = tmp5_0*w15;
2721                                        const double tmp39_1 = A_10_2*w7;
2722                                        const double tmp40_1 = tmp11_0*w17;
2723                                        const double tmp41_1 = A_01_3*w2;
2724                                        const double tmp42_1 = tmp12_0*w16;
2725                                        const double tmp43_1 = tmp17_0*w5;
2726                                        const double tmp44_1 = tmp7_0*w7;
2727                                        const double tmp45_1 = tmp6_0*w4;
2728                                        const double tmp46_1 = A_01_3*w8;
2729                                        const double tmp47_1 = A_10_2*w4;
2730                                        const double tmp48_1 = A_10_1*w7;
2731                                        const double tmp49_1 = tmp12_0*w17;
2732                                        const double tmp50_1 = A_01_0*w2;
2733                                        const double tmp51_1 = tmp11_0*w16;
2734                                        const double tmp52_1 = tmp18_0*w8;
2735                                        const double tmp53_1 = tmp19_0*w2;
2736                                        const double tmp54_1 = tmp13_0*w8;
2737                                        const double tmp55_1 = tmp15_0*w2;
2738                                        const double tmp56_1 = tmp19_0*w8;
2739                                        const double tmp57_1 = tmp18_0*w2;
2740                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
2741                                        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;
2742                                        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;
2743                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
2744                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
2745                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
2746                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
2747                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
2748                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2749                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
2750                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
2751                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
2752                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
2753                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
2754                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
2755                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2756                                    }
2757                                }
2758                            } else { // constant data
2759                                for (index_t k=0; k<numEq; k++) {
2760                                    for (index_t m=0; m<numComp; m++) {
2761                                        const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
2762                                        const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
2763                                        const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
2764                                        const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
2765                                        const double tmp0_0 = A_01 + A_10;
2766                                        const double tmp0_1 = A_00*w18;
2767                                        const double tmp1_1 = A_01*w19;
2768                                        const double tmp2_1 = A_10*w20;
2769                                        const double tmp3_1 = A_11*w21;
2770                                        const double tmp4_1 = A_00*w22;
2771                                        const double tmp5_1 = tmp0_0*w19;
2772                                        const double tmp6_1 = A_11*w23;
2773                                        const double tmp7_1 = A_11*w25;
2774                                        const double tmp8_1 = A_00*w24;
2775                                        const double tmp9_1 = tmp0_0*w20;
2776                                        const double tmp10_1 = A_01*w20;
2777                                        const double tmp11_1 = A_11*w27;
2778                                        const double tmp12_1 = A_00*w26;
2779                                        const double tmp13_1 = A_10*w19;
2780                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2781