/[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 3781 by caltinay, Mon Jan 23 22:18:49 2012 UTC trunk/ripley/src/Rectangle.cpp revision 4013 by caltinay, Thu Oct 4 03:13:27 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2012 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
17  extern "C" {  extern "C" {
18  #include "paso/SystemMatrix.h"  #include <paso/SystemMatrix.h>
19  }  }
20    
21    #ifdef USE_NETCDF
22    #include <netcdfcpp.h>
23    #endif
24    
25  #if USE_SILO  #if USE_SILO
26  #include <silo.h>  #include <silo.h>
27  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 32  namespace ripley { Line 38  namespace ripley {
38  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
39                       double y1, int d0, int d1) :                       double y1, int d0, int d1) :
40      RipleyDomain(2),      RipleyDomain(2),
     m_gNE0(n0),  
     m_gNE1(n1),  
41      m_x0(x0),      m_x0(x0),
42      m_y0(y0),      m_y0(y0),
43      m_l0(x1-x0),      m_l0(x1-x0),
44      m_l1(y1-y0),      m_l1(y1-y0)
     m_NX(d0),  
     m_NY(d1)  
45  {  {
46        // ignore subdivision parameters for serial run
47        if (m_mpiInfo->size == 1) {
48            d0=1;
49            d1=1;
50        }
51    
52        bool warn=false;
53        // if number of subdivisions is non-positive, try to subdivide by the same
54        // ratio as the number of elements
55        if (d0<=0 && d1<=0) {
56            warn=true;
57            d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
58            d1=m_mpiInfo->size/d0;
59            if (d0*d1 != m_mpiInfo->size) {
60                // ratios not the same so subdivide side with more elements only
61                if (n0>n1) {
62                    d0=0;
63                    d1=1;
64                } else {
65                    d0=1;
66                    d1=0;
67                }
68            }
69        }
70        if (d0<=0) {
71            // d1 is preset, determine d0 - throw further down if result is no good
72            d0=m_mpiInfo->size/d1;
73        } else if (d1<=0) {
74            // d0 is preset, determine d1 - throw further down if result is no good
75            d1=m_mpiInfo->size/d0;
76        }
77    
78        m_NX=d0;
79        m_NY=d1;
80    
81      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
82      // among number of ranks      // among number of ranks
83      if (m_NX*m_NY != m_mpiInfo->size)      if (m_NX*m_NY != m_mpiInfo->size)
84          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
85    
86      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)      if (warn) {
87          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
88                << d1 << "). This may not be optimal!" << endl;
89        }
90    
91        if ((n0+1)%m_NX > 0) {
92            double Dx=m_l0/n0;
93            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
94            m_l0=Dx*n0;
95            cout << "Warning: Adjusted number of elements and length. N0="
96                << n0 << ", l0=" << m_l0 << endl;
97        }
98        if ((n1+1)%m_NY > 0) {
99            double Dy=m_l1/n1;
100            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
101            m_l1=Dy*n1;
102            cout << "Warning: Adjusted number of elements and length. N1="
103                << n1 << ", l1=" << m_l1 << endl;
104        }
105    
106        m_gNE0=n0;
107        m_gNE1=n1;
108    
109      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
110          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
# Line 83  Rectangle::Rectangle(int n0, int n1, dou Line 140  Rectangle::Rectangle(int n0, int n1, dou
140    
141  Rectangle::~Rectangle()  Rectangle::~Rectangle()
142  {  {
143        Paso_SystemMatrixPattern_free(m_pattern);
144        Paso_Connector_free(m_connector);
145  }  }
146    
147  string Rectangle::getDescription() const  string Rectangle::getDescription() const
# Line 104  bool Rectangle::operator==(const Abstrac Line 163  bool Rectangle::operator==(const Abstrac
163      return false;      return false;
164  }  }
165    
166    void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
167                const vector<int>& first, const vector<int>& numValues) const
168    {
169    #ifdef USE_NETCDF
170        // check destination function space
171        int myN0, myN1;
172        if (out.getFunctionSpace().getTypeCode() == Nodes) {
173            myN0 = m_N0;
174            myN1 = m_N1;
175        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
176                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
177            myN0 = m_NE0;
178            myN1 = m_NE1;
179        } else
180            throw RipleyException("readNcGrid(): invalid function space for output data object");
181    
182        if (first.size() != 2)
183            throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
184    
185        if (numValues.size() != 2)
186            throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
187    
188        // check file existence and size
189        NcFile f(filename.c_str(), NcFile::ReadOnly);
190        if (!f.is_valid())
191            throw RipleyException("readNcGrid(): cannot open file");
192    
193        NcVar* var = f.get_var(varname.c_str());
194        if (!var)
195            throw RipleyException("readNcGrid(): invalid variable");
196    
197        // TODO: rank>0 data support
198        const int numComp = out.getDataPointSize();
199        if (numComp > 1)
200            throw RipleyException("readNcGrid(): only scalar data supported");
201    
202        const int dims = var->num_dims();
203        const long *edges = var->edges();
204    
205        // is this a slice of the data object (dims!=2)?
206        // note the expected ordering of edges (as in numpy: y,x)
207        if ( (dims==2 && (numValues[1] > edges[0] || numValues[0] > edges[1]))
208                || (dims==1 && numValues[1]>1) ) {
209            throw RipleyException("readNcGrid(): not enough data in file");
210        }
211    
212        // check if this rank contributes anything
213        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
214                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1)
215            return;
216    
217        // now determine how much this rank has to write
218    
219        // first coordinates in data object to write to
220        const int first0 = max(0, first[0]-m_offset0);
221        const int first1 = max(0, first[1]-m_offset1);
222        // indices to first value in file
223        const int idx0 = max(0, m_offset0-first[0]);
224        const int idx1 = max(0, m_offset1-first[1]);
225        // number of values to write
226        const int num0 = min(numValues[0]-idx0, myN0-first0);
227        const int num1 = min(numValues[1]-idx1, myN1-first1);
228    
229        vector<double> values(num0*num1);
230        if (dims==2) {
231            var->set_cur(idx1, idx0);
232            var->get(&values[0], num1, num0);
233        } else {
234            var->set_cur(idx0);
235            var->get(&values[0], num0);
236        }
237    
238        const int dpp = out.getNumDataPointsPerSample();
239        out.requireWrite();
240    
241        for (index_t y=0; y<num1; y++) {
242    #pragma omp parallel for
243            for (index_t x=0; x<num0; x++) {
244                const int dataIndex = (first1+y)*myN0+first0+x;
245                const int srcIndex=y*num0+x;
246                double* dest = out.getSampleDataRW(dataIndex);
247                for (index_t q=0; q<dpp; q++) {
248                    *dest++ = values[srcIndex];
249                }
250            }
251        }
252    #else
253        throw RipleyException("readNcGrid(): not compiled with netCDF support");
254    #endif
255    }
256    
257    void Rectangle::readBinaryGrid(escript::Data& out, string filename,
258                const vector<int>& first, const vector<int>& numValues) const
259    {
260        // check destination function space
261        int myN0, myN1;
262        if (out.getFunctionSpace().getTypeCode() == Nodes) {
263            myN0 = m_N0;
264            myN1 = m_N1;
265        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
266                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
267            myN0 = m_NE0;
268            myN1 = m_NE1;
269        } else
270            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
271    
272        // check file existence and size
273        ifstream f(filename.c_str(), ifstream::binary);
274        if (f.fail()) {
275            throw RipleyException("readBinaryGrid(): cannot open file");
276        }
277        f.seekg(0, ios::end);
278        const int numComp = out.getDataPointSize();
279        const int filesize = f.tellg();
280        const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(float);
281        if (filesize < reqsize) {
282            f.close();
283            throw RipleyException("readBinaryGrid(): not enough data in file");
284        }
285    
286        // check if this rank contributes anything
287        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
288                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1) {
289            f.close();
290            return;
291        }
292    
293        // now determine how much this rank has to write
294    
295        // first coordinates in data object to write to
296        const int first0 = max(0, first[0]-m_offset0);
297        const int first1 = max(0, first[1]-m_offset1);
298        // indices to first value in file
299        const int idx0 = max(0, m_offset0-first[0]);
300        const int idx1 = max(0, m_offset1-first[1]);
301        // number of values to write
302        const int num0 = min(numValues[0]-idx0, myN0-first0);
303        const int num1 = min(numValues[1]-idx1, myN1-first1);
304    
305        out.requireWrite();
306        vector<float> values(num0*numComp);
307        const int dpp = out.getNumDataPointsPerSample();
308    
309        for (index_t y=0; y<num1; y++) {
310            const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);
311            f.seekg(fileofs*sizeof(float));
312            f.read((char*)&values[0], num0*numComp*sizeof(float));
313            for (index_t x=0; x<num0; x++) {
314                double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0);
315                for (index_t c=0; c<numComp; c++) {
316                    for (index_t q=0; q<dpp; q++) {
317                        *dest++ = static_cast<double>(values[x*numComp+c]);
318                    }
319                }
320            }
321        }
322    
323        f.close();
324    }
325    
326  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
327  {  {
328  #if USE_SILO  #if USE_SILO
# Line 112  void Rectangle::dump(const string& fileN Line 331  void Rectangle::dump(const string& fileN
331          fn+=".silo";          fn+=".silo";
332      }      }
333    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
334      int driver=DB_HDF5;          int driver=DB_HDF5;    
335      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
336        const char* blockDirFmt = "/block%04d";
337    
338  #ifdef ESYS_MPI  #ifdef ESYS_MPI
339      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
340        const int NUM_SILO_FILES = 1;
341  #endif  #endif
342    
343      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 244  void Rectangle::dump(const string& fileN Line 463  void Rectangle::dump(const string& fileN
463      }      }
464    
465  #else // USE_SILO  #else // USE_SILO
466      throw RipleyException("dump(): no Silo support");      throw RipleyException("dump: no Silo support");
467  #endif  #endif
468  }  }
469    
# Line 268  const int* Rectangle::borrowSampleRefere Line 487  const int* Rectangle::borrowSampleRefere
487      }      }
488    
489      stringstream msg;      stringstream msg;
490      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
491      throw RipleyException(msg.str());      throw RipleyException(msg.str());
492  }  }
493    
# Line 318  bool Rectangle::ownSample(int fsType, in Line 536  bool Rectangle::ownSample(int fsType, in
536      }      }
537    
538      stringstream msg;      stringstream msg;
539      msg << "ownSample() not implemented for "      msg << "ownSample: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
540      throw RipleyException(msg.str());      throw RipleyException(msg.str());
541  }  }
542    
# Line 420  void Rectangle::setToNormal(escript::Dat Line 637  void Rectangle::setToNormal(escript::Dat
637    
638      } else {      } else {
639          stringstream msg;          stringstream msg;
640          msg << "setToNormal() not implemented for "          msg << "setToNormal: invalid function space type "
641              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());              << out.getFunctionSpace().getTypeCode();
642          throw RipleyException(msg.str());          throw RipleyException(msg.str());
643      }      }
644  }  }
# Line 434  void Rectangle::setToSize(escript::Data& Line 651  void Rectangle::setToSize(escript::Data&
651          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
652          const double hSize=getFirstCoordAndSpacing(0).second;          const double hSize=getFirstCoordAndSpacing(0).second;
653          const double vSize=getFirstCoordAndSpacing(1).second;          const double vSize=getFirstCoordAndSpacing(1).second;
654          const double size=min(hSize,vSize);          const double size=sqrt(hSize*hSize+vSize*vSize);
655  #pragma omp parallel for  #pragma omp parallel for
656          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
657              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 483  void Rectangle::setToSize(escript::Data& Line 700  void Rectangle::setToSize(escript::Data&
700    
701      } else {      } else {
702          stringstream msg;          stringstream msg;
703          msg << "setToSize() not implemented for "          msg << "setToSize: invalid function space type "
704              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());              << out.getFunctionSpace().getTypeCode();
705          throw RipleyException(msg.str());          throw RipleyException(msg.str());
706      }      }
707  }  }
# Line 565  pair<double,double> Rectangle::getFirstC Line 782  pair<double,double> Rectangle::getFirstC
782      } else if (dim==1) {      } else if (dim==1) {
783          return pair<double,double>(m_y0+(m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);          return pair<double,double>(m_y0+(m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
784      }      }
785      throw RipleyException("getFirstCoordAndSpacing(): invalid argument");      throw RipleyException("getFirstCoordAndSpacing: invalid argument");
786  }  }
787    
788  //protected  //protected
# Line 632  void Rectangle::assembleGradient(escript Line 849  void Rectangle::assembleGradient(escript
849    
850      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
851          out.requireWrite();          out.requireWrite();
852  #pragma omp parallel for  #pragma omp parallel
853          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
854              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
855                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
856                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
857                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
858                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
859                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
860                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
861                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
862                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
863                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
864                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
865                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
866                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      for (index_t i=0; i < numComp; ++i) {
867                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
868                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
869                  } /* end of component loop i */                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
870              } /* end of k0 loop */                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
871          } /* end of k1 loop */                          o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
872                            o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
873                            o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
874                            o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
875                        } // end of component loop i
876                    } // end of k0 loop
877                } // end of k1 loop
878            } // end of parallel section
879      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
880          out.requireWrite();          out.requireWrite();
881  #pragma omp parallel for  #pragma omp parallel
882          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
883              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
884                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
885                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
886                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
887                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
888                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
889                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
890                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
891                      o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
892                  } /* end of component loop i */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
893              } /* end of k0 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
894          } /* end of k1 loop */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
895                        for (index_t i=0; i < numComp; ++i) {
896                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
897                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
898                        } // end of component loop i
899                    } // end of k0 loop
900                } // end of k1 loop
901            } // end of parallel section
902      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
903          out.requireWrite();          out.requireWrite();
904  #pragma omp parallel  #pragma omp parallel
905          {          {
906                vector<double> f_00(numComp);
907                vector<double> f_01(numComp);
908                vector<double> f_10(numComp);
909                vector<double> f_11(numComp);
910              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
911  #pragma omp for nowait  #pragma omp for nowait
912                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
913                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
914                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
915                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
916                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
917                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
918                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
919                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
920                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
921                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
922                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
923                      } /* end of component loop i */                      } // end of component loop i
924                  } /* end of k1 loop */                  } // end of k1 loop
925              } /* end of face 0 */              } // end of face 0
926              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
927  #pragma omp for nowait  #pragma omp for nowait
928                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
929                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
930                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
931                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
932                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
933                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
934                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
935                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
936                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
937                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
938                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
939                      } /* end of component loop i */                      } // end of component loop i
940                  } /* end of k1 loop */                  } // end of k1 loop
941              } /* end of face 1 */              } // end of face 1
942              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
943  #pragma omp for nowait  #pragma omp for nowait
944                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
945                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
946                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
947                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
948                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
949                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
950                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
951                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
952                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
953                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
954                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
955                      } /* end of component loop i */                      } // end of component loop i
956                  } /* end of k0 loop */                  } // end of k0 loop
957              } /* end of face 2 */              } // end of face 2
958              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
959  #pragma omp for nowait  #pragma omp for nowait
960                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
961                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
962                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
963                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
964                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
965                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
966                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
967                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
968                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
969                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
970                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
971                      } /* end of component loop i */                      } // end of component loop i
972                  } /* end of k0 loop */                  } // end of k0 loop
973              } /* end of face 3 */              } // end of face 3
974          } // end of parallel section          } // end of parallel section
975    
976      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
977          out.requireWrite();          out.requireWrite();
978  #pragma omp parallel  #pragma omp parallel
979          {          {
980                vector<double> f_00(numComp);
981                vector<double> f_01(numComp);
982                vector<double> f_10(numComp);
983                vector<double> f_11(numComp);
984              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
985  #pragma omp for nowait  #pragma omp for nowait
986                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
987                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
988                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
989                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
990                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
991                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
992                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
993                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
994                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
995                      } /* end of component loop i */                      } // end of component loop i
996                  } /* end of k1 loop */                  } // end of k1 loop
997              } /* end of face 0 */              } // end of face 0
998              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
999  #pragma omp for nowait  #pragma omp for nowait
1000                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1001                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
1002                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
1003                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
1004                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
1005                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1006                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1007                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
1008                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
1009                      } /* end of component loop i */                      } // end of component loop i
1010                  } /* end of k1 loop */                  } // end of k1 loop
1011              } /* end of face 1 */              } // end of face 1
1012              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1013  #pragma omp for nowait  #pragma omp for nowait
1014                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1015                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
1016                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
1017                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
1018                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
1019                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1020                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1021                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
1022                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
1023                      } /* end of component loop i */                      } // end of component loop i
1024                  } /* end of k0 loop */                  } // end of k0 loop
1025              } /* end of face 2 */              } // end of face 2
1026              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1027  #pragma omp for nowait  #pragma omp for nowait
1028                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1029                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
1030                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
1031                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
1032                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
1033                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1034                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1035                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
1036                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
1037                      } /* end of component loop i */                      } // end of component loop i
1038                  } /* end of k0 loop */                  } // end of k0 loop
1039              } /* end of face 3 */              } // end of face 3
1040          } // end of parallel section          } // end of parallel section
1041      }      }
1042  }  }
# Line 809  void Rectangle::assembleIntegrate(vector Line 1049  void Rectangle::assembleIntegrate(vector
1049      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
1050      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1051      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1052      if (arg.getFunctionSpace().getTypeCode() == Elements) {      const int fs=arg.getFunctionSpace().getTypeCode();
1053          const double w = h0*h1/4.;      if (fs == Elements && arg.actsExpanded()) {
1054  #pragma omp parallel  #pragma omp parallel
1055          {          {
1056              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1057                const double w = h0*h1/4.;
1058  #pragma omp for nowait  #pragma omp for nowait
1059              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1060                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
# Line 824  void Rectangle::assembleIntegrate(vector Line 1065  void Rectangle::assembleIntegrate(vector
1065                          const double f2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
1066                          const double f3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
1067                          int_local[i]+=(f0+f1+f2+f3)*w;                          int_local[i]+=(f0+f1+f2+f3)*w;
1068                      }  /* end of component loop i */                      }  // end of component loop i
1069                  } /* end of k0 loop */                  } // end of k0 loop
1070              } /* end of k1 loop */              } // end of k1 loop
   
1071  #pragma omp critical  #pragma omp critical
1072              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1073                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1074          } // end of parallel section          } // end of parallel section
1075      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
1076        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1077          const double w = h0*h1;          const double w = h0*h1;
1078  #pragma omp parallel  #pragma omp parallel
1079          {          {
# Line 843  void Rectangle::assembleIntegrate(vector Line 1084  void Rectangle::assembleIntegrate(vector
1084                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
1085                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1086                          int_local[i]+=f[i]*w;                          int_local[i]+=f[i]*w;
1087                      }  /* end of component loop i */                      }
1088                  } /* end of k0 loop */                  }
1089              } /* end of k1 loop */              }
   
1090  #pragma omp critical  #pragma omp critical
1091              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1092                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1093          } // end of parallel section          } // end of parallel section
1094      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
1095          const double w0 = h0/2.;      } else if (fs == FaceElements && arg.actsExpanded()) {
         const double w1 = h1/2.;  
1096  #pragma omp parallel  #pragma omp parallel
1097          {          {
1098              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1099                const double w0 = h0/2.;
1100                const double w1 = h1/2.;
1101              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1102  #pragma omp for nowait  #pragma omp for nowait
1103                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
# Line 865  void Rectangle::assembleIntegrate(vector Line 1106  void Rectangle::assembleIntegrate(vector
1106                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1107                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1108                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1109                      }  /* end of component loop i */                      }  // end of component loop i
1110                  } /* end of k1 loop */                  } // end of k1 loop
1111              }              }
1112    
1113              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
# Line 877  void Rectangle::assembleIntegrate(vector Line 1118  void Rectangle::assembleIntegrate(vector
1118                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1119                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1120                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1121                      }  /* end of component loop i */                      }  // end of component loop i
1122                  } /* end of k1 loop */                  } // end of k1 loop
1123              }              }
1124    
1125              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
# Line 889  void Rectangle::assembleIntegrate(vector Line 1130  void Rectangle::assembleIntegrate(vector
1130                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1131                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1132                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1133                      }  /* end of component loop i */                      }  // end of component loop i
1134                  } /* end of k0 loop */                  } // end of k0 loop
1135              }              }
1136    
1137              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
# Line 901  void Rectangle::assembleIntegrate(vector Line 1142  void Rectangle::assembleIntegrate(vector
1142                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1143                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1144                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1145                      }  /* end of component loop i */                      }  // end of component loop i
1146                  } /* end of k0 loop */                  } // end of k0 loop
1147              }              }
   
1148  #pragma omp critical  #pragma omp critical
1149              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1150                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1151          } // end of parallel section          } // end of parallel section
1152      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
1153        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1154  #pragma omp parallel  #pragma omp parallel
1155          {          {
1156              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
# Line 919  void Rectangle::assembleIntegrate(vector Line 1160  void Rectangle::assembleIntegrate(vector
1160                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1161                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1162                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
1163                      }  /* end of component loop i */                      }
1164                  } /* end of k1 loop */                  }
1165              }              }
1166    
1167              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
# Line 929  void Rectangle::assembleIntegrate(vector Line 1170  void Rectangle::assembleIntegrate(vector
1170                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1171                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1172                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
1173                      }  /* end of component loop i */                      }
1174                  } /* end of k1 loop */                  }
1175              }              }
1176    
1177              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
# Line 939  void Rectangle::assembleIntegrate(vector Line 1180  void Rectangle::assembleIntegrate(vector
1180                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1181                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1182                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
1183                      }  /* end of component loop i */                      }
1184                  } /* end of k0 loop */                  }
1185              }              }
1186    
1187              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
# Line 949  void Rectangle::assembleIntegrate(vector Line 1190  void Rectangle::assembleIntegrate(vector
1190                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1191                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1192                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
1193                      }  /* end of component loop i */                      }
1194                  } /* end of k0 loop */                  }
1195              }              }
1196    
1197  #pragma omp critical  #pragma omp critical
1198              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1199                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1200          } // end of parallel section          } // end of parallel section
1201      }      } // function space selector
1202  }  }
1203    
1204  //protected  //protected
# Line 1027  void Rectangle::dofToNodes(escript::Data Line 1268  void Rectangle::dofToNodes(escript::Data
1268                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1269          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1270      }      }
1271        Paso_Coupler_free(coupler);
1272  }  }
1273    
1274  //private  //private
# Line 1309  void Rectangle::interpolateNodesOnElemen Line 1551  void Rectangle::interpolateNodesOnElemen
1551      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1552      if (reduced) {      if (reduced) {
1553          out.requireWrite();          out.requireWrite();
1554          const double c0 = .25;          const double c0 = 0.25;
1555  #pragma omp parallel for  #pragma omp parallel
1556          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1557              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1558                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_01(numComp);
1559                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1560                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              vector<double> f_11(numComp);
1561                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1562                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1563                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1564                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
1565                  } /* end of component loop i */                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
1566              } /* end of k0 loop */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
1567          } /* end of k1 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
1568                        double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1569                        for (index_t i=0; i < numComp; ++i) {
1570                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1571                        } /* end of component loop i */
1572                    } /* end of k0 loop */
1573                } /* end of k1 loop */
1574            } /* end of parallel section */
1575      } else {      } else {
1576          out.requireWrite();          out.requireWrite();
1577          const double c0 = .16666666666666666667;          const double c0 = 0.16666666666666666667;
1578          const double c1 = .044658198738520451079;          const double c1 = 0.044658198738520451079;
1579          const double c2 = .62200846792814621559;          const double c2 = 0.62200846792814621559;
1580  #pragma omp parallel for  #pragma omp parallel
1581          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1582              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1583                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
1584                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1585                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
1586                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1587                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1588                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1589                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
1590                      o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
1591                      o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
1592                      o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
1593                  } /* end of component loop i */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1594              } /* end of k0 loop */                      for (index_t i=0; i < numComp; ++i) {
1595          } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1596                            o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1597                            o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1598                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1599                        } /* end of component loop i */
1600                    } /* end of k0 loop */
1601                } /* end of k1 loop */
1602            } /* end of parallel section */
1603      }      }
1604  }  }
1605    
# Line 1354  void Rectangle::interpolateNodesOnFaces( Line 1610  void Rectangle::interpolateNodesOnFaces(
1610      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1611      if (reduced) {      if (reduced) {
1612          out.requireWrite();          out.requireWrite();
1613          const double c0 = .5;          const double c0 = 0.5;
1614  #pragma omp parallel  #pragma omp parallel
1615          {          {
1616                vector<double> f_00(numComp);
1617                vector<double> f_01(numComp);
1618                vector<double> f_10(numComp);
1619                vector<double> f_11(numComp);
1620              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1621  #pragma omp for nowait  #pragma omp for nowait
1622                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1623                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
1624                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
1625                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1626                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1627                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
# Line 1371  void Rectangle::interpolateNodesOnFaces( Line 1631  void Rectangle::interpolateNodesOnFaces(
1631              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1632  #pragma omp for nowait  #pragma omp for nowait
1633                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1634                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
1635                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
1636                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1637                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1638                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
# Line 1382  void Rectangle::interpolateNodesOnFaces( Line 1642  void Rectangle::interpolateNodesOnFaces(
1642              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1643  #pragma omp for nowait  #pragma omp for nowait
1644                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1645                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
1646                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
1647                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1648                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1649                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
# Line 1393  void Rectangle::interpolateNodesOnFaces( Line 1653  void Rectangle::interpolateNodesOnFaces(
1653              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1654  #pragma omp for nowait  #pragma omp for nowait
1655                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1656                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
1657                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
1658                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1659                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1660                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1661                      } /* end of component loop i */                      } /* end of component loop i */
1662                  } /* end of k0 loop */                  } /* end of k0 loop */
1663              } /* end of face 3 */              } /* end of face 3 */
1664          } // end of parallel section          } /* end of parallel section */
1665      } else {      } else {
1666          out.requireWrite();          out.requireWrite();
1667          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1668          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1669  #pragma omp parallel  #pragma omp parallel
1670          {          {
1671                vector<double> f_00(numComp);
1672                vector<double> f_01(numComp);
1673                vector<double> f_10(numComp);
1674                vector<double> f_11(numComp);
1675              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1676  #pragma omp for nowait  #pragma omp for nowait
1677                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1678                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
1679                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
1680                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1681                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1682                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_01[i] + c1*f_00[i];
1683                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_01[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_00[i] + c1*f_01[i];
1684                      } /* end of component loop i */                      } /* end of component loop i */
1685                  } /* end of k1 loop */                  } /* end of k1 loop */
1686              } /* end of face 0 */              } /* end of face 0 */
1687              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1688  #pragma omp for nowait  #pragma omp for nowait
1689                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1690                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
1691                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
1692                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1693                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1694                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = c1*f_10[i] + c0*f_11[i];
1695                          o[INDEX2(i,numComp,1)] = f_10[i]*c0 + f_11[i]*c1;                          o[INDEX2(i,numComp,1)] = c1*f_11[i] + c0*f_10[i];
1696                      } /* end of component loop i */                      } /* end of component loop i */
1697                  } /* end of k1 loop */                  } /* end of k1 loop */
1698              } /* end of face 1 */              } /* end of face 1 */
1699              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1700  #pragma omp for nowait  #pragma omp for nowait
1701                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1702                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
1703                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
1704                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1705                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1706                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_10[i] + c1*f_00[i];
1707                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_10[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_00[i] + c1*f_10[i];
1708                      } /* end of component loop i */                      } /* end of component loop i */
1709                  } /* end of k0 loop */                  } /* end of k0 loop */
1710              } /* end of face 2 */              } /* end of face 2 */
1711              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1712  #pragma omp for nowait  #pragma omp for nowait
1713                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1714                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
1715                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
1716                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1717                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1718                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_11[i] + c1*f_01[i];
1719                          o[INDEX2(i,numComp,1)] = f_01[i]*c0 + f_11[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_01[i] + c1*f_11[i];
1720                      } /* end of component loop i */                      } /* end of component loop i */
1721                  } /* end of k0 loop */                  } /* end of k0 loop */
1722              } /* end of face 3 */              } /* end of face 3 */
1723          } // end of parallel section          } /* end of parallel section */
1724      }      }
1725  }  }
1726    
# Line 2729  void Rectangle::assemblePDESystem(Paso_S Line 2993  void Rectangle::assemblePDESystem(Paso_S
2993                                      const double tmp78_1 = B_0_0*w34;                                      const double tmp78_1 = B_0_0*w34;
2994                                      const double tmp12_1 = tmp0_0*w36;                                      const double tmp12_1 = tmp0_0*w36;
2995                                      const double tmp75_1 = B_1_0*w32;                                      const double tmp75_1 = B_1_0*w32;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
2996                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
2997                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2998                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2999                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
3000                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
3001                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
3002                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
3003                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
3004                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
3005                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
3006                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
3007                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
3008                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
3009                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
3010                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
3011                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
3012                                  }                                  }
3013                              }                              }
3014                          } else { // constant data                          } else { // constant data
# Line 2760  void Rectangle::assemblePDESystem(Paso_S Line 3024  void Rectangle::assemblePDESystem(Paso_S
3024                                      const double tmp2_1 = B_0*w46;                                      const double tmp2_1 = B_0*w46;
3025                                      const double tmp7_1 = B_0*w51;                                      const double tmp7_1 = B_0*w51;
3026                                      const double tmp3_1 = B_0*w47;                                      const double tmp3_1 = B_0*w47;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
3027                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;
3028                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3029                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;
3030                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
3031                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3032                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;
3033                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;
3034                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;
3035                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;
3036                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;  
3037                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;
3038                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
3039                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;
3040                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;
3041                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;
3042                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;
3043                                  }                                  }
3044                              }                              }
3045                          }                          }
# Line 3353  void Rectangle::assemblePDEBoundarySingl Line 3617  void Rectangle::assemblePDEBoundarySingl
3617  {  {
3618      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
3619      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
     /* GENERATOR SNIP_PDEBC_SINGLE_PRE TOP */  
3620      const double w0 = 0.31100423396407310779*h1;      const double w0 = 0.31100423396407310779*h1;
3621      const double w1 = 0.022329099369260225539*h1;      const double w1 = 0.022329099369260225539*h1;
3622      const double w10 = 0.022329099369260225539*h0;      const double w10 = 0.022329099369260225539*h0;
# Line 3370  void Rectangle::assemblePDEBoundarySingl Line 3633  void Rectangle::assemblePDEBoundarySingl
3633      const double w7 = 0.5*h1;      const double w7 = 0.5*h1;
3634      const double w8 = 0.083333333333333333333*h0;      const double w8 = 0.083333333333333333333*h0;
3635      const double w9 = 0.31100423396407310779*h0;      const double w9 = 0.31100423396407310779*h0;
     /* GENERATOR SNIP_PDEBC_SINGLE_PRE BOTTOM */  
3636      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3637      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3638      rhs.requireWrite();      rhs.requireWrite();
# Line 3378  void Rectangle::assemblePDEBoundarySingl Line 3640  void Rectangle::assemblePDEBoundarySingl
3640      {      {
3641          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3642              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3643  #pragma omp for nowait  #pragma omp for
3644                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3645                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3646                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3647                      const index_t e = k1;                      const index_t e = k1;
                     /* GENERATOR SNIP_PDEBC_SINGLE_0 TOP */  
3648                      ///////////////                      ///////////////
3649                      // process d //                      // process d //
3650                      ///////////////                      ///////////////
# Line 3399  void Rectangle::assemblePDEBoundarySingl Line 3660  void Rectangle::assemblePDEBoundarySingl
3660                              const double tmp3_1 = d_1*w0;                              const double tmp3_1 = d_1*w0;
3661                              const double tmp2_1 = tmp0_0*w2;                              const double tmp2_1 = tmp0_0*w2;
3662                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
                             EM_S[INDEX2(0,2,4)]+=tmp2_1;  
3663                              EM_S[INDEX2(2,0,4)]+=tmp2_1;                              EM_S[INDEX2(2,0,4)]+=tmp2_1;
3664                                EM_S[INDEX2(0,2,4)]+=tmp2_1;
3665                              EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;
3666                          } else { /* constant data */                          } else { /* constant data */
3667                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3668                              const double tmp1_1 = d_0*w4;                              const double tmp1_1 = d_0*w4;
3669                              const double tmp0_1 = d_0*w3;                              const double tmp0_1 = d_0*w3;
3670                              EM_S[INDEX2(0,0,4)]+=tmp0_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1;
                             EM_S[INDEX2(0,2,4)]+=tmp1_1;  
3671                              EM_S[INDEX2(2,0,4)]+=tmp1_1;                              EM_S[INDEX2(2,0,4)]+=tmp1_1;
3672                                EM_S[INDEX2(0,2,4)]+=tmp1_1;
3673                              EM_S[INDEX2(2,2,4)]+=tmp0_1;                              EM_S[INDEX2(2,2,4)]+=tmp0_1;
3674                          }                          }
3675                      }                      }
# Line 3433  void Rectangle::assemblePDEBoundarySingl Line 3694  void Rectangle::assemblePDEBoundarySingl
3694                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
3695                          }                          }
3696                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_0 BOTTOM */  
3697                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_N0*k1;
3698                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3699                  }                  }
# Line 3442  void Rectangle::assemblePDEBoundarySingl Line 3702  void Rectangle::assemblePDEBoundarySingl
3702    
3703          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
3704              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3705  #pragma omp for nowait  #pragma omp for
3706                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3707                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3708                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3709                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
                     /* GENERATOR SNIP_PDEBC_SINGLE_1 TOP */  
3710                      ///////////////                      ///////////////
3711                      // process d //                      // process d //
3712                      ///////////////                      ///////////////
# Line 3462  void Rectangle::assemblePDEBoundarySingl Line 3721  void Rectangle::assemblePDEBoundarySingl
3721                              const double tmp3_1 = d_0*w0;                              const double tmp3_1 = d_0*w0;
3722                              const double tmp0_1 = d_1*w0;                              const double tmp0_1 = d_1*w0;
3723                              const double tmp2_1 = tmp0_0*w2;                              const double tmp2_1 = tmp0_0*w2;
3724                              EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;
3725                              EM_S[INDEX2(3,1,4)]+=tmp2_1;                              EM_S[INDEX2(3,1,4)]+=tmp2_1;
3726                              EM_S[INDEX2(1,3,4)]+=tmp2_1;                              EM_S[INDEX2(1,3,4)]+=tmp2_1;
3727                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp1_1;
3728                          } else { /* constant data */                          } else { /* constant data */
3729                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3730                              const double tmp1_1 = d_0*w4;                              const double tmp1_1 = d_0*w4;
3731                              const double tmp0_1 = d_0*w3;                              const double tmp0_1 = d_0*w3;
3732                              EM_S[INDEX2(3,3,4)]+=tmp0_1;                              EM_S[INDEX2(1,1,4)]+=tmp0_1;
3733                              EM_S[INDEX2(3,1,4)]+=tmp1_1;                              EM_S[INDEX2(3,1,4)]+=tmp1_1;
3734                              EM_S[INDEX2(1,3,4)]+=tmp1_1;                              EM_S[INDEX2(1,3,4)]+=tmp1_1;
3735                              EM_S[INDEX2(1,1,4)]+=tmp0_1;                              EM_S[INDEX2(3,3,4)]+=tmp0_1;
3736                          }                          }
3737                      }                      }
3738                      ///////////////                      ///////////////
# Line 3497  void Rectangle::assemblePDEBoundarySingl Line 3756  void Rectangle::assemblePDEBoundarySingl
3756                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
3757                          }                          }
3758                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_1 BOTTOM */  
3759                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_N0*(k1+1)-2;
3760                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3761                  }                  }
# Line 3506  void Rectangle::assemblePDEBoundarySingl Line 3764  void Rectangle::assemblePDEBoundarySingl
3764    
3765          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
3766              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3767  #pragma omp for nowait  #pragma omp for
3768                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3769                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3770                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3771                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
                     /* GENERATOR SNIP_PDEBC_SINGLE_2 TOP */  
3772                      ///////////////                      ///////////////
3773                      // process d //                      // process d //
3774                      ///////////////                      ///////////////
# Line 3526  void Rectangle::assemblePDEBoundarySingl Line 3783  void Rectangle::assemblePDEBoundarySingl
3783                              const double tmp0_1 = tmp0_0*w8;                              const double tmp0_1 = tmp0_0*w8;
3784                              const double tmp1_1 = d_1*w10;                              const double tmp1_1 = d_1*w10;
3785                              const double tmp3_1 = d_0*w10;                              const double tmp3_1 = d_0*w10;
                             EM_S[INDEX2(0,1,4)]+=tmp0_1;  
3786                              EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;                              EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;
3787                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
3788                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
3789                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;
3790                          } else { /* constant data */                          } else { /* constant data */
3791                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3792                              const double tmp0_1 = d_0*w11;                              const double tmp0_1 = d_0*w11;
3793                              const double tmp1_1 = d_0*w12;                              const double tmp1_1 = d_0*w12;
                             EM_S[INDEX2(0,1,4)]+=tmp0_1;  
3794                              EM_S[INDEX2(0,0,4)]+=tmp1_1;                              EM_S[INDEX2(0,0,4)]+=tmp1_1;
3795                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
3796                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
3797                              EM_S[INDEX2(1,1,4)]+=tmp1_1;                              EM_S[INDEX2(1,1,4)]+=tmp1_1;
3798                          }                          }
3799                      }                      }
# Line 3561  void Rectangle::assemblePDEBoundarySingl Line 3818  void Rectangle::assemblePDEBoundarySingl
3818                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
3819                          }                          }
3820                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_2 BOTTOM */  
3821                      const index_t firstNode=k0;                      const index_t firstNode=k0;
3822                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3823                  }                  }
# Line 3570  void Rectangle::assemblePDEBoundarySingl Line 3826  void Rectangle::assemblePDEBoundarySingl
3826    
3827          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
3828              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3829  #pragma omp for nowait  #pragma omp for
3830                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3831                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
3832                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3833                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
                     /* GENERATOR SNIP_PDEBC_SINGLE_3 TOP */  
3834                      ///////////////                      ///////////////
3835                      // process d //                      // process d //
3836                      ///////////////                      ///////////////
# Line 3590  void Rectangle::assemblePDEBoundarySingl Line 3845  void Rectangle::assemblePDEBoundarySingl
3845                              const double tmp0_1 = tmp0_0*w8;                              const double tmp0_1 = tmp0_0*w8;
3846                              const double tmp3_1 = d_1*w10;                              const double tmp3_1 = d_1*w10;
3847                              const double tmp1_1 = d_0*w10;                              const double tmp1_1 = d_0*w10;
3848                                EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;
3849                              EM_S[INDEX2(3,2,4)]+=tmp0_1;                              EM_S[INDEX2(3,2,4)]+=tmp0_1;
                             EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;  
3850                              EM_S[INDEX2(2,3,4)]+=tmp0_1;                              EM_S[INDEX2(2,3,4)]+=tmp0_1;
3851                              EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
3852                          } else { /* constant data */                          } else { /* constant data */
3853                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3854                              const double tmp0_1 = d_0*w11;                              const double tmp0_1 = d_0*w11;
3855                              const double tmp1_1 = d_0*w12;                              const double tmp1_1 = d_0*w12;
3856                                EM_S[INDEX2(2,2,4)]+=tmp1_1;
3857                              EM_S[INDEX2(3,2,4)]+=tmp0_1;                              EM_S[INDEX2(3,2,4)]+=tmp0_1;
                             EM_S[INDEX2(3,3,4)]+=tmp1_1;  
3858                              EM_S[INDEX2(2,3,4)]+=tmp0_1;                              EM_S[INDEX2(2,3,4)]+=tmp0_1;
3859                              EM_S[INDEX2(2,2,4)]+=tmp1_1;                              EM_S[INDEX2(3,3,4)]+=tmp1_1;
3860                          }                          }
3861                      }                      }
3862                      ///////////////                      ///////////////
# Line 3625  void Rectangle::assemblePDEBoundarySingl Line 3880  void Rectangle::assemblePDEBoundarySingl
3880                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
3881                          }                          }
3882                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_3 BOTTOM */  
3883                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_N0*(m_N1-2)+k0;
3884                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3885                  }                  }
# Line 3640  void Rectangle::assemblePDEBoundarySingl Line 3894  void Rectangle::assemblePDEBoundarySingl
3894  {  {
3895      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
3896      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE TOP */  
3897      const double w0 = 0.25*h1;      const double w0 = 0.25*h1;
3898      const double w1 = 0.5*h1;      const double w1 = 0.5*h1;
3899      const double w2 = 0.25*h0;      const double w2 = 0.25*h0;
3900      const double w3 = 0.5*h0;      const double w3 = 0.5*h0;
     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE BOTTOM */  
3901      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3902      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3903      rhs.requireWrite();      rhs.requireWrite();
# Line 3653  void Rectangle::assemblePDEBoundarySingl Line 3905  void Rectangle::assemblePDEBoundarySingl
3905      {      {
3906          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3907              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3908  #pragma omp for nowait  #pragma omp for
3909                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3910                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3911                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3912                      const index_t e = k1;                      const index_t e = k1;
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 TOP */  
3913                      ///////////////                      ///////////////
3914                      // process d //                      // process d //
3915                      ///////////////                      ///////////////
# Line 3667  void Rectangle::assemblePDEBoundarySingl Line 3918  void Rectangle::assemblePDEBoundarySingl
3918                          const double d_0 = d_p[0];                          const double d_0 = d_p[0];
3919                          const double tmp0_1 = d_0*w0;                          const double tmp0_1 = d_0*w0;
3920                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1;
                         EM_S[INDEX2(0,2,4)]+=tmp0_1;  
3921                          EM_S[INDEX2(2,0,4)]+=tmp0_1;                          EM_S[INDEX2(2,0,4)]+=tmp0_1;
3922                            EM_S[INDEX2(0,2,4)]+=tmp0_1;
3923                          EM_S[INDEX2(2,2,4)]+=tmp0_1;                          EM_S[INDEX2(2,2,4)]+=tmp0_1;
3924                      }                      }
3925                      ///////////////                      ///////////////
# Line 3676  void Rectangle::assemblePDEBoundarySingl Line 3927  void Rectangle::assemblePDEBoundarySingl
3927                      ///////////////                      ///////////////
3928                      if (add_EM_F) {                      if (add_EM_F) {
3929                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3930                          const double y_0 = y_p[0];                          const double tmp0_1 = w1*y_p[0];
                         const double tmp0_1 = w1*y_0;  
3931                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0_1;
3932                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
3933                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 BOTTOM */  
3934                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_N0*k1;
3935                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3936                  }                  }
# Line 3690  void Rectangle::assemblePDEBoundarySingl Line 3939  void Rectangle::assemblePDEBoundarySingl
3939    
3940          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
3941              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3942  #pragma omp for nowait  #pragma omp for
3943                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3944                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3945                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3946                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 TOP */  
3947                      ///////////////                      ///////////////
3948                      // process d //                      // process d //
3949                      ///////////////                      ///////////////
# Line 3703  void Rectangle::assemblePDEBoundarySingl Line 3951  void Rectangle::assemblePDEBoundarySingl
3951                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3952                          const double d_0 = d_p[0];                          const double d_0 = d_p[0];
3953                          const double tmp0_1 = d_0*w0;                          const double tmp0_1 = d_0*w0;
3954                          EM_S[INDEX2(3,3,4)]+=tmp0_1;                          EM_S[INDEX2(1,1,4)]+=tmp0_1;
3955                          EM_S[INDEX2(3,1,4)]+=tmp0_1;                          EM_S[INDEX2(3,1,4)]+=tmp0_1;
3956                          EM_S[INDEX2(1,3,4)]+=tmp0_1;                          EM_S[INDEX2(1,3,4)]+=tmp0_1;
3957                          EM_S[INDEX2(1,1,4)]+=tmp0_1;                          EM_S[INDEX2(3,3,4)]+=tmp0_1;
3958                      }                      }
3959                      ///////////////                      ///////////////
3960                      // process y //                      // process y //
3961                      ///////////////                      ///////////////
3962                      if (add_EM_F) {                      if (add_EM_F) {
3963                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3964                          const double y_0 = y_p[0];                          const double tmp0_1 = w1*y_p[0];
                         const double tmp0_1 = w1*y_0;  
3965                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0_1;
3966                          EM_F[3]+=tmp0_1;                          EM_F[3]+=tmp0_1;
3967                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 BOTTOM */  
3968                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_N0*(k1+1)-2;
3969                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3970                  }                  }
# Line 3727  void Rectangle::assemblePDEBoundarySingl Line 3973  void Rectangle::assemblePDEBoundarySingl
3973    
3974          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
3975              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3976  #pragma omp for nowait  #pragma omp for
3977                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3978                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3979                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3980                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 TOP */  
3981                      ///////////////                      ///////////////
3982                      // process d //                      // process d //
3983                      ///////////////                      ///////////////
3984                      if (add_EM_S) {                      if (add_EM_S) {
3985                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3986                          const double d_0 = d_p[0];                          const double tmp0_1 = d_p[0]*w2;
                         const double tmp0_1 = d_0*w2;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1;  
3987                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1;
3988                          EM_S[INDEX2(1,0,4)]+=tmp0_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1;
3989                            EM_S[INDEX2(0,1,4)]+=tmp0_1;
3990                          EM_S[INDEX2(1,1,4)]+=tmp0_1;                          EM_S[INDEX2(1,1,4)]+=tmp0_1;
3991                      }                      }
3992                      ///////////////                      ///////////////
# Line 3750  void Rectangle::assemblePDEBoundarySingl Line 3994  void Rectangle::assemblePDEBoundarySingl
3994                      ///////////////                      ///////////////
3995                      if (add_EM_F) {                      if (add_EM_F) {
3996                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3997                          const double y_0 = y_p[0];                          const double tmp0_1 = w3*y_p[0];
                         const double tmp0_1 = w3*y_0;  
3998                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0_1;
3999                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0_1;
4000                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 BOTTOM */  
4001                      const index_t firstNode=k0;                      const index_t firstNode=k0;
4002                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
4003                  }                  }
# Line 3764  void Rectangle::assemblePDEBoundarySingl Line 4006  void Rectangle::assemblePDEBoundarySingl
4006    
4007          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4008              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4009  #pragma omp for nowait  #pragma omp for
4010                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4011                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
4012                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
4013                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 TOP */  
4014                      ///////////////                      ///////////////
4015                      // process d //                      // process d //
4016                      ///////////////                      ///////////////
4017                      if (add_EM_S) {                      if (add_EM_S) {
4018                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
4019                          const double d_0 = d_p[0];                          const double tmp0_1 = d_p[0]*w2;
4020                          const double tmp0_1 = d_0*w2;                          EM_S[INDEX2(2,2,4)]+=tmp0_1;
4021                          EM_S[INDEX2(3,2,4)]+=tmp0_1;                          EM_S[INDEX2(3,2,4)]+=tmp0_1;
                         EM_S[INDEX2(3,3,4)]+=tmp0_1;  
4022                          EM_S[INDEX2(2,3,4)]+=tmp0_1;                          EM_S[INDEX2(2,3,4)]+=tmp0_1;
4023                          EM_S[INDEX2(2,2,4)]+=tmp0_1;                          EM_S[INDEX2(3,3,4)]+=tmp0_1;
4024                      }                      }
4025                      ///////////////                      ///////////////
4026                      // process y //                      // process y //
4027                      ///////////////                      ///////////////
4028                      if (add_EM_F) {                      if (add_EM_F) {
4029                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4030                          const double y_0 = y_p[0];                          const double tmp0_1 = w3*y_p[0];
                         const double tmp0_1 = w3*y_0;  
4031                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
4032                          EM_F[3]+=tmp0_1;                          EM_F[3]+=tmp0_1;
4033                      }                      }
                     /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 BOTTOM */  
4034                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_N0*(m_N1-2)+k0;
4035                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
4036                  }                  }
# Line 3814  void Rectangle::assemblePDEBoundarySyste Line 4052  void Rectangle::assemblePDEBoundarySyste
4052          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
4053          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
4054      }      }
     /* GENERATOR SNIP_PDEBC_SYSTEM_PRE TOP */  
4055      const double w0 = 0.31100423396407310779*h1;      const double w0 = 0.31100423396407310779*h1;
4056      const double w1 = 0.022329099369260225539*h1;      const double w1 = 0.022329099369260225539*h1;
4057      const double w10 = 0.022329099369260225539*h0;      const double w10 = 0.022329099369260225539*h0;
# Line 3831  void Rectangle::assemblePDEBoundarySyste Line 4068  void Rectangle::assemblePDEBoundarySyste
4068      const double w7 = 0.5*h1;      const double w7 = 0.5*h1;
4069      const double w8 = 0.083333333333333333333*h0;      const double w8 = 0.083333333333333333333*h0;
4070      const double w9 = 0.31100423396407310779*h0;      const double w9 = 0.31100423396407310779*h0;
     /* GENERATOR SNIP_PDEBC_SYSTEM_PRE BOTTOM */  
4071      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
4072      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
4073      rhs.requireWrite();      rhs.requireWrite();
# Line 3839  void Rectangle::assemblePDEBoundarySyste Line 4075  void Rectangle::assemblePDEBoundarySyste
4075      {      {
4076          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
4077              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
4078  #pragma omp for nowait  #pragma omp for
4079                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4080                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4081                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4082                      const index_t e = k1;                      const index_t e = k1;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_0 TOP */  
4083                      ///////////////                      ///////////////
4084                      // process d //                      // process d //
4085                      ///////////////                      ///////////////
# Line 3856  void Rectangle::assemblePDEBoundarySyste Line 4091  void Rectangle::assemblePDEBoundarySyste
4091                                      const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
4092                                      const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
4093                                      const double tmp0_0 = d_0 + d_1;                                      const double tmp0_0 = d_0 + d_1;
                                     const double tmp1_1 = d_1*w1;  
                                     const double tmp4_1 = d_0*w1;  
4094                                      const double tmp0_1 = d_0*w0;                                      const double tmp0_1 = d_0*w0;
4095                                      const double tmp3_1 = d_1*w0;                                      const double tmp1_1 = d_1*w1;
4096                                      const double tmp2_1 = tmp0_0*w2;                                      const double tmp2_1 = tmp0_0*w2;
4097                                        const double tmp3_1 = d_1*w0;
4098                                        const double tmp4_1 = d_0*w1;
4099                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1;  
4100                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1;
4101                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1;
4102                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
4103                                  }                                  }
4104                               }                               }
# Line 3871  void Rectangle::assemblePDEBoundarySyste Line 4106  void Rectangle::assemblePDEBoundarySyste
4106                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
4107                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
4108                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
                                     const double tmp1_1 = d_0*w4;  
4109                                      const double tmp0_1 = d_0*w3;                                      const double tmp0_1 = d_0*w3;
4110                                        const double tmp1_1 = d_0*w4;
4111                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp1_1;  
4112                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1;
4113                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp1_1;
4114                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
4115                                  }                                  }
4116                              }                              }
# Line 3906  void Rectangle::assemblePDEBoundarySyste Line 4141  void Rectangle::assemblePDEBoundarySyste
4141                              }                              }
4142                          }                          }
4143                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_0 BOTTOM */  
4144                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_N0*k1;
4145                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4146                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 3916  void Rectangle::assemblePDEBoundarySyste Line 4150  void Rectangle::assemblePDEBoundarySyste
4150    
4151          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
4152              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4153  #pragma omp for nowait  #pragma omp for
4154                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4155                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4156                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4157                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_1 TOP */  
4158                      ///////////////                      ///////////////
4159                      // process d //                      // process d //
4160                      ///////////////                      ///////////////
# Line 3938  void Rectangle::assemblePDEBoundarySyste Line 4171  void Rectangle::assemblePDEBoundarySyste
4171                                      const double tmp3_1 = d_0*w0;                                      const double tmp3_1 = d_0*w0;
4172                                      const double tmp0_1 = d_1*w0;                                      const double tmp0_1 = d_1*w0;
4173                                      const double tmp2_1 = tmp0_0*w2;                                      const double tmp2_1 = tmp0_0*w2;
4174                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
4175                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1;
4176                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1;
4177                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
4178                                  }                                  }
4179                               }                               }
4180                          } else { /* constant data */                          } else { /* constant data */
# Line 3950  void Rectangle::assemblePDEBoundarySyste Line 4183  void Rectangle::assemblePDEBoundarySyste
4183                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
4184                                      const double tmp1_1 = d_0*w4;                                      const double tmp1_1 = d_0*w4;
4185                                      const double tmp0_1 = d_0*w3;                                      const double tmp0_1 = d_0*w3;
4186                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
4187                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp1_1;
4188                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1;
4189                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
4190                                  }                                  }
4191                              }                              }
4192                          }                          }
# Line 3983  void Rectangle::assemblePDEBoundarySyste Line 4216  void Rectangle::assemblePDEBoundarySyste
4216                              }                              }
4217                          }                          }
4218                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_1 BOTTOM */  
4219                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_N0*(k1+1)-2;
4220                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4221                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 3993  void Rectangle::assemblePDEBoundarySyste Line 4225  void Rectangle::assemblePDEBoundarySyste
4225    
4226          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
4227              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4228  #pragma omp for nowait  #pragma omp for
4229                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4230                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4231                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4232                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_2 TOP */  
4233                      ///////////////                      ///////////////
4234                      // process d //                      // process d //
4235                      ///////////////                      ///////////////
# Line 4015  void Rectangle::assemblePDEBoundarySyste Line 4246  void Rectangle::assemblePDEBoundarySyste
4246                                      const double tmp0_1 = tmp0_0*w8;                                      const double tmp0_1 = tmp0_0*w8;
4247                                      const double tmp1_1 = d_1*w10;                                      const double tmp1_1 = d_1*w10;
4248                                      const double tmp3_1 = d_0*w10;                                      const double tmp3_1 = d_0*w10;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
4249                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
4250                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
4251                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
4252                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
4253                                  }                                  }
4254                               }                               }
# Line 4027  void Rectangle::assemblePDEBoundarySyste Line 4258  void Rectangle::assemblePDEBoundarySyste
4258                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
4259                                      const double tmp0_1 = d_0*w11;                                      const double tmp0_1 = d_0*w11;
4260                                      const double tmp1_1 = d_0*w12;                                      const double tmp1_1 = d_0*w12;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
4261                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1;
4262                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
4263                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
4264                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1;
4265                                  }                                  }
4266                              }                              }
# Line 4060  void Rectangle::assemblePDEBoundarySyste Line 4291  void Rectangle::assemblePDEBoundarySyste
4291                              }                              }
4292                          }                          }
4293                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_2 BOTTOM */  
4294                      const index_t firstNode=k0;                      const index_t firstNode=k0;
4295                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4296                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 4070  void Rectangle::assemblePDEBoundarySyste Line 4300  void Rectangle::assemblePDEBoundarySyste
4300    
4301          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4302              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4303  #pragma omp for nowait  #pragma omp for
4304                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4305                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4306                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4307                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_3 TOP */  
4308                      ///////////////                      ///////////////
4309                      // process d //                      // process d //
4310                      ///////////////                      ///////////////
# Line 4092  void Rectangle::assemblePDEBoundarySyste Line 4321  void Rectangle::assemblePDEBoundarySyste
4321                                      const double tmp0_1 = tmp0_0*w8;                                      const double tmp0_1 = tmp0_0*w8;
4322                                      const double tmp3_1 = d_1*w10;                                      const double tmp3_1 = d_1*w10;
4323                                      const double tmp1_1 = d_0*w10;                                      const double tmp1_1 = d_0*w10;
4324                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
4325                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
4326                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
4327                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
4328                                  }                                  }
4329                               }                               }
4330                          } else { /* constant data */                          } else { /* constant data */
# Line 4104  void Rectangle::assemblePDEBoundarySyste Line 4333  void Rectangle::assemblePDEBoundarySyste
4333                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
4334                                      const double tmp0_1 = d_0*w11;                                      const double tmp0_1 = d_0*w11;
4335                                      const double tmp1_1 = d_0*w12;                                      const double tmp1_1 = d_0*w12;
4336                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp1_1;
4337                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1;  
4338                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
4339                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1;
4340                                  }                                  }
4341                              }                              }
4342                          }                          }
# Line 4137  void Rectangle::assemblePDEBoundarySyste Line 4366  void Rectangle::assemblePDEBoundarySyste
4366                              }                              }
4367                          }                          }
4368                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_3 BOTTOM */  
4369                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_N0*(m_N1-2)+k0;
4370                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4371                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 4160  void Rectangle::assemblePDEBoundarySyste Line 4388  void Rectangle::assemblePDEBoundarySyste
4388          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
4389          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
4390      }      }
     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_PRE TOP */  
4391      const double w0 = 0.25*h1;      const double w0 = 0.25*h1;
4392      const double w1 = 0.5*h1;      const double w1 = 0.5*h1;
4393      const double w2 = 0.25*h0;      const double w2 = 0.25*h0;
4394      const double w3 = 0.5*h0;      const double w3 = 0.5*h0;
     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_PRE BOTTOM */  
4395      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
4396      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
4397      rhs.requireWrite();      rhs.requireWrite();
# Line 4173  void Rectangle::assemblePDEBoundarySyste Line 4399  void Rectangle::assemblePDEBoundarySyste
4399      {      {
4400          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
4401              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
4402  #pragma omp for nowait  #pragma omp for
4403                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4404                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4405                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4406                      const index_t e = k1;                      const index_t e = k1;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_0 TOP */  
4407                      ///////////////                      ///////////////
4408                      // process d //                      // process d //
4409                      ///////////////                      ///////////////
# Line 4189  void Rectangle::assemblePDEBoundarySyste Line 4414  void Rectangle::assemblePDEBoundarySyste
4414                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double d_0 = d_p[INDEX2(k, m, numEq)];
4415                                  const double tmp0_1 = d_0*w0;                                  const double tmp0_1 = d_0*w0;
4416                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
                                 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;  
4417                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;
4418                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;
4419                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
4420                              }                              }
4421                          }                          }
# Line 4201  void Rectangle::assemblePDEBoundarySyste Line 4426  void Rectangle::assemblePDEBoundarySyste
4426                      if (add_EM_F) {                      if (add_EM_F) {
4427                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4428                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4429                              const double y_0 = y_p[k];                              const double tmp0_1 = w1*y_p[k];
                             const double tmp0_1 = w1*y_0;  
4430                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
4431                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
4432                          }                          }
4433                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_0 BOTTOM */  
4434                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_N0*k1;
4435                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4436                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
4437                  }                  }
                 // ADD EM_F  and EM_S  
4438              } // end colouring              } // end colouring
4439          }          }
4440    
4441          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
4442              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4443  #pragma omp for nowait  #pragma omp for
4444                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4445                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4446                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4447                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_1 TOP */  
4448                      ///////////////                      ///////////////
4449                      // process d //                      // process d //
4450                      ///////////////                      ///////////////
# Line 4233  void Rectangle::assemblePDEBoundarySyste Line 4454  void Rectangle::assemblePDEBoundarySyste
4454                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
4455                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double d_0 = d_p[INDEX2(k, m, numEq)];
4456                                  const double tmp0_1 = d_0*w0;                                  const double tmp0_1 = d_0*w0;
4457                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
4458                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;
4459                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;
4460                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
4461                              }                              }
4462                          }                          }
4463                      }                      }
# Line 4246  void Rectangle::assemblePDEBoundarySyste Line 4467  void Rectangle::assemblePDEBoundarySyste
4467                      if (add_EM_F) {                      if (add_EM_F) {
4468                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4469                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4470                              const double y_0 = y_p[k];                              const double tmp0_1 = w1*y_p[k];
                             const double tmp0_1 = w1*y_0;  
4471                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
4472                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
4473                          }                          }
4474                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_1 BOTTOM */  
4475                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_N0*(k1+1)-2;
4476                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4477                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 4262  void Rectangle::assemblePDEBoundarySyste Line 4481  void Rectangle::assemblePDEBoundarySyste
4481    
4482          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
4483              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4484  #pragma omp for nowait  #pragma omp for
4485                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4486                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4487                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4488                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_2 TOP */  
4489                      ///////////////                      ///////////////
4490                      // process d //                      // process d //
4491                      ///////////////                      ///////////////
# Line 4277  void Rectangle::assemblePDEBoundarySyste Line 4495  void Rectangle::assemblePDEBoundarySyste
4495                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
4496                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double d_0 = d_p[INDEX2(k, m, numEq)];
4497                                  const double tmp0_1 = d_0*w2;                                  const double tmp0_1 = d_0*w2;
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
4498                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
4499                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
4500                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
4501                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
4502                              }                              }
4503                          }                          }
# Line 4290  void Rectangle::assemblePDEBoundarySyste Line 4508  void Rectangle::assemblePDEBoundarySyste
4508                      if (add_EM_F) {                      if (add_EM_F) {
4509                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4510                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4511                              const double y_0 = y_p[k];                              const double tmp0_1 = w3*y_p[k];
                             const double tmp0_1 = w3*y_0;  
4512                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
4513                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
4514                          }                          }
4515                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_2 BOTTOM */  
4516                      const index_t firstNode=k0;                      const index_t firstNode=k0;
4517                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4518                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
# Line 4306  void Rectangle::assemblePDEBoundarySyste Line 4522  void Rectangle::assemblePDEBoundarySyste
4522    
4523          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4524              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4525  #pragma omp for nowait  #pragma omp for
4526                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4527                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4528                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4529                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_3 TOP */  
4530                      ///////////////                      ///////////////
4531                      // process d //                      // process d //
4532                      ///////////////                      ///////////////
# Line 4321  void Rectangle::assemblePDEBoundarySyste Line 4536  void Rectangle::assemblePDEBoundarySyste
4536                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
4537                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double d_0 = d_p[INDEX2(k, m, numEq)];
4538                                  const double tmp0_1 = d_0*w2;                                  const double tmp0_1 = d_0*w2;
4539                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
4540                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;  
4541                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
4542                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
4543                              }                              }
4544                          }                          }
4545                      }                      }
# Line 4334  void Rectangle::assemblePDEBoundarySyste Line 4549  void Rectangle::assemblePDEBoundarySyste
4549                      if (add_EM_F) {                      if (add_EM_F) {
4550                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4551                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4552                              const double y_0 = y_p[k];                              const double tmp0_1 = w3*y_p[k];
                             const double tmp0_1 = w3*y_0;  
4553                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
4554                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
4555                          }                          }
4556                      }                      }
                     /* GENERATOR SNIP_PDEBC_SYSTEM_REDUCED_3 BOTTOM */  
4557                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_N0*(m_N1-2)+k0;
4558                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
4559                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);

Legend:
Removed from v.3781  
changed lines
  Added in v.4013

  ViewVC Help
Powered by ViewVC 1.1.26