/[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

revision 3793 by gross, Wed Feb 1 07:39:43 2012 UTC revision 4368 by jfenwick, Thu Apr 18 02:26:35 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2012 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
 extern "C" {  
17  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
18  }  #include <esysUtils/esysFileWriter.h>
19    
20    #ifdef USE_NETCDF
21    #include <netcdfcpp.h>
22    #endif
23    
24  #if USE_SILO  #if USE_SILO
25  #include <silo.h>  #include <silo.h>
# Line 26  extern "C" { Line 31  extern "C" {
31  #include <iomanip>  #include <iomanip>
32    
33  using namespace std;  using namespace std;
34    using esysUtils::FileWriter;
35    
36  namespace ripley {  namespace ripley {
37    
38  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
39                       double y1, int d0, int d1) :                       double y1, int d0, int d1) :
40      RipleyDomain(2),      RipleyDomain(2)
     m_gNE0(n0),  
     m_gNE1(n1),  
     m_x0(x0),  
     m_y0(y0),  
     m_l0(x1-x0),  
     m_l1(y1-y0),  
     m_NX(d0),  
     m_NY(d1)  
41  {  {
42        // ignore subdivision parameters for serial run
43        if (m_mpiInfo->size == 1) {
44            d0=1;
45            d1=1;
46        }
47    
48        bool warn=false;
49        // if number of subdivisions is non-positive, try to subdivide by the same
50        // ratio as the number of elements
51        if (d0<=0 && d1<=0) {
52            warn=true;
53            d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
54            d1=m_mpiInfo->size/d0;
55            if (d0*d1 != m_mpiInfo->size) {
56                // ratios not the same so subdivide side with more elements only
57                if (n0>n1) {
58                    d0=0;
59                    d1=1;
60                } else {
61                    d0=1;
62                    d1=0;
63                }
64            }
65        }
66        if (d0<=0) {
67            // d1 is preset, determine d0 - throw further down if result is no good
68            d0=m_mpiInfo->size/d1;
69        } else if (d1<=0) {
70            // d0 is preset, determine d1 - throw further down if result is no good
71            d1=m_mpiInfo->size/d0;
72        }
73    
74      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
75      // among number of ranks      // among number of ranks
76      if (m_NX*m_NY != m_mpiInfo->size)      if (d0*d1 != m_mpiInfo->size)
77          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
78    
79      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)      if (warn) {
80          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
81                << d1 << "). This may not be optimal!" << endl;
82        }
83    
84        double l0 = x1-x0;
85        double l1 = y1-y0;
86        m_dx[0] = l0/n0;
87        m_dx[1] = l1/n1;
88    
89        if ((n0+1)%d0 > 0) {
90            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
91            l0=m_dx[0]*n0;
92            cout << "Warning: Adjusted number of elements and length. N0="
93                << n0 << ", l0=" << l0 << endl;
94        }
95        if ((n1+1)%d1 > 0) {
96            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
97            l1=m_dx[1]*n1;
98            cout << "Warning: Adjusted number of elements and length. N1="
99                << n1 << ", l1=" << l1 << endl;
100        }
101    
102      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))      if ((d0 > 1 && (n0+1)/d0<2) || (d1 > 1 && (n1+1)/d1<2))
103          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
104    
105        m_gNE[0] = n0;
106        m_gNE[1] = n1;
107        m_origin[0] = x0;
108        m_origin[1] = y0;
109        m_length[0] = l0;
110        m_length[1] = l1;
111        m_NX[0] = d0;
112        m_NX[1] = d1;
113    
114      // local number of elements (with and without overlap)      // local number of elements (with and without overlap)
115      m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);      m_NE[0] = m_ownNE[0] = (d0>1 ? (n0+1)/d0 : n0);
116      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)      if (m_mpiInfo->rank%d0>0 && m_mpiInfo->rank%d0<d0-1)
117          m_NE0++;          m_NE[0]++;
118      else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)      else if (d0>1 && m_mpiInfo->rank%d0==d0-1)
119          m_ownNE0--;          m_ownNE[0]--;
120    
121      m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);      m_NE[1] = m_ownNE[1] = (d1>1 ? (n1+1)/d1 : n1);
122      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)      if (m_mpiInfo->rank/d0>0 && m_mpiInfo->rank/d0<d1-1)
123          m_NE1++;          m_NE[1]++;
124      else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)      else if (d1>1 && m_mpiInfo->rank/d0==d1-1)
125          m_ownNE1--;          m_ownNE[1]--;
126    
127      // local number of nodes      // local number of nodes
128      m_N0 = m_NE0+1;      m_NN[0] = m_NE[0]+1;
129      m_N1 = m_NE1+1;      m_NN[1] = m_NE[1]+1;
130    
131      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
132      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);      m_offset[0] = (n0+1)/d0*(m_mpiInfo->rank%d0);
133      if (m_offset0 > 0)      if (m_offset[0] > 0)
134          m_offset0--;          m_offset[0]--;
135      m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);      m_offset[1] = (n1+1)/d1*(m_mpiInfo->rank/d0);
136      if (m_offset1 > 0)      if (m_offset[1] > 0)
137          m_offset1--;          m_offset[1]--;
138    
139      populateSampleIds();      populateSampleIds();
140      createPattern();      createPattern();
# Line 97  bool Rectangle::operator==(const Abstrac Line 156  bool Rectangle::operator==(const Abstrac
156      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
157      if (o) {      if (o) {
158          return (RipleyDomain::operator==(other) &&          return (RipleyDomain::operator==(other) &&
159                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1                  m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1]
160                  && m_x0==o->m_x0 && m_y0==o->m_y0                  && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1]
161                  && m_l0==o->m_l0 && m_l1==o->m_l1                  && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1]
162                  && m_NX==o->m_NX && m_NY==o->m_NY);                  && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1]);
163      }      }
164    
165      return false;      return false;
166  }  }
167    
168    void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
169                const vector<int>& first, const vector<int>& numValues,
170                const vector<int>& multiplier) const
171    {
172    #ifdef USE_NETCDF
173        // check destination function space
174        int myN0, myN1;
175        if (out.getFunctionSpace().getTypeCode() == Nodes) {
176            myN0 = m_NN[0];
177            myN1 = m_NN[1];
178        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
179                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
180            myN0 = m_NE[0];
181            myN1 = m_NE[1];
182        } else
183            throw RipleyException("readNcGrid(): invalid function space for output data object");
184    
185        if (first.size() != 2)
186            throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
187    
188        if (numValues.size() != 2)
189            throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
190    
191        if (multiplier.size() != 2)
192            throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");
193        for (size_t i=0; i<multiplier.size(); i++)
194            if (multiplier[i]<1)
195                throw RipleyException("readNcGrid(): all multipliers must be positive");
196    
197        // check file existence and size
198        NcFile f(filename.c_str(), NcFile::ReadOnly);
199        if (!f.is_valid())
200            throw RipleyException("readNcGrid(): cannot open file");
201    
202        NcVar* var = f.get_var(varname.c_str());
203        if (!var)
204            throw RipleyException("readNcGrid(): invalid variable");
205    
206        // TODO: rank>0 data support
207        const int numComp = out.getDataPointSize();
208        if (numComp > 1)
209            throw RipleyException("readNcGrid(): only scalar data supported");
210    
211        const int dims = var->num_dims();
212        const long *edges = var->edges();
213    
214        // is this a slice of the data object (dims!=2)?
215        // note the expected ordering of edges (as in numpy: y,x)
216        if ( (dims==2 && (numValues[1] > edges[0] || numValues[0] > edges[1]))
217                || (dims==1 && numValues[1]>1) ) {
218            throw RipleyException("readNcGrid(): not enough data in file");
219        }
220    
221        // check if this rank contributes anything
222        if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset[0] ||
223                first[1] >= m_offset[1]+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset[1])
224            return;
225    
226        // now determine how much this rank has to write
227    
228        // first coordinates in data object to write to
229        const int first0 = max(0, first[0]-m_offset[0]);
230        const int first1 = max(0, first[1]-m_offset[1]);
231        // indices to first value in file
232        const int idx0 = max(0, m_offset[0]-first[0]);
233        const int idx1 = max(0, m_offset[1]-first[1]);
234        // number of values to read
235        const int num0 = min(numValues[0]-idx0, myN0-first0);
236        const int num1 = min(numValues[1]-idx1, myN1-first1);
237    
238        vector<double> values(num0*num1);
239        if (dims==2) {
240            var->set_cur(idx1, idx0);
241            var->get(&values[0], num1, num0);
242        } else {
243            var->set_cur(idx0);
244            var->get(&values[0], num0);
245        }
246    
247        const int dpp = out.getNumDataPointsPerSample();
248        out.requireWrite();
249    
250        for (index_t y=0; y<num1; y++) {
251    #pragma omp parallel for
252            for (index_t x=0; x<num0; x++) {
253                const int baseIndex = first0+x*multiplier[0]
254                                      +(first1+y*multiplier[1])*myN0;
255                const int srcIndex = y*num0+x;
256                if (!isnan(values[srcIndex])) {
257                    for (index_t m1=0; m1<multiplier[1]; m1++) {
258                        for (index_t m0=0; m0<multiplier[0]; m0++) {
259                            const int dataIndex = baseIndex+m0+m1*myN0;
260                            double* dest = out.getSampleDataRW(dataIndex);
261                            for (index_t q=0; q<dpp; q++) {
262                                *dest++ = values[srcIndex];
263                            }
264                        }
265                    }
266                }
267            }
268        }
269    #else
270        throw RipleyException("readNcGrid(): not compiled with netCDF support");
271    #endif
272    }
273    
274    void Rectangle::readBinaryGrid(escript::Data& out, string filename,
275                                   const vector<int>& first,
276                                   const vector<int>& numValues,
277                                   const vector<int>& multiplier) const
278    {
279        // check destination function space
280        int myN0, myN1;
281        if (out.getFunctionSpace().getTypeCode() == Nodes) {
282            myN0 = m_NN[0];
283            myN1 = m_NN[1];
284        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
285                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
286            myN0 = m_NE[0];
287            myN1 = m_NE[1];
288        } else
289            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
290    
291        // check file existence and size
292        ifstream f(filename.c_str(), ifstream::binary);
293        if (f.fail()) {
294            throw RipleyException("readBinaryGrid(): cannot open file");
295        }
296        f.seekg(0, ios::end);
297        const int numComp = out.getDataPointSize();
298        const int filesize = f.tellg();
299        const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(float);
300        if (filesize < reqsize) {
301            f.close();
302            throw RipleyException("readBinaryGrid(): not enough data in file");
303        }
304    
305        // check if this rank contributes anything
306        if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0] <= m_offset[0] ||
307                first[1] >= m_offset[1]+myN1 || first[1]+numValues[1] <= m_offset[1]) {
308            f.close();
309            return;
310        }
311    
312        // now determine how much this rank has to write
313    
314        // first coordinates in data object to write to
315        const int first0 = max(0, first[0]-m_offset[0]);
316        const int first1 = max(0, first[1]-m_offset[1]);
317        // indices to first value in file
318        const int idx0 = max(0, m_offset[0]-first[0]);
319        const int idx1 = max(0, m_offset[1]-first[1]);
320        // number of values to read
321        const int num0 = min(numValues[0]-idx0, myN0-first0);
322        const int num1 = min(numValues[1]-idx1, myN1-first1);
323    
324        out.requireWrite();
325        vector<float> values(num0*numComp);
326        const int dpp = out.getNumDataPointsPerSample();
327    
328        for (index_t y=0; y<num1; y++) {
329            const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);
330            f.seekg(fileofs*sizeof(float));
331            f.read((char*)&values[0], num0*numComp*sizeof(float));
332            for (index_t x=0; x<num0; x++) {
333                const int baseIndex = first0+x*multiplier[0]
334                                        +(first1+y*multiplier[1])*myN0;
335                for (index_t m1=0; m1<multiplier[1]; m1++) {
336                    for (index_t m0=0; m0<multiplier[0]; m0++) {
337                        const int dataIndex = baseIndex+m0+m1*myN0;
338                        double* dest = out.getSampleDataRW(dataIndex);
339                        for (index_t c=0; c<numComp; c++) {
340                            if (!std::isnan(values[x*numComp+c])) {
341                                for (index_t q=0; q<dpp; q++) {
342                                    *dest++ = static_cast<double>(values[x*numComp+c]);
343                                }
344                            }
345                        }
346                    }
347                }
348            }
349        }
350    
351        f.close();
352    }
353    
354    void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,
355                                    int byteOrder, int dataType) const
356    {
357        // the mapping is not universally correct but should work on our
358        // supported platforms
359        switch (dataType) {
360            case DATATYPE_INT32:
361                writeBinaryGridImpl<int>(in, filename, byteOrder);
362                break;
363            case DATATYPE_FLOAT32:
364                writeBinaryGridImpl<float>(in, filename, byteOrder);
365                break;
366            case DATATYPE_FLOAT64:
367                writeBinaryGridImpl<double>(in, filename, byteOrder);
368                break;
369            default:
370                throw RipleyException("writeBinaryGrid(): invalid or unsupported datatype");
371        }
372    }
373    
374    template<typename ValueType>
375    void Rectangle::writeBinaryGridImpl(const escript::Data& in,
376                                        const string& filename, int byteOrder) const
377    {
378        // check function space and determine number of points
379        int myN0, myN1;
380        int totalN0, totalN1;
381        if (in.getFunctionSpace().getTypeCode() == Nodes) {
382            myN0 = m_NN[0];
383            myN1 = m_NN[1];
384            totalN0 = m_gNE[0]+1;
385            totalN1 = m_gNE[1]+1;
386        } else if (in.getFunctionSpace().getTypeCode() == Elements ||
387                    in.getFunctionSpace().getTypeCode() == ReducedElements) {
388            myN0 = m_NE[0];
389            myN1 = m_NE[1];
390            totalN0 = m_gNE[0];
391            totalN1 = m_gNE[1];
392        } else
393            throw RipleyException("writeBinaryGrid(): invalid function space of data object");
394    
395        const int numComp = in.getDataPointSize();
396        const int dpp = in.getNumDataPointsPerSample();
397    
398        if (numComp > 1 || dpp > 1)
399            throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
400    
401        escript::Data* _in = const_cast<escript::Data*>(&in);
402        const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;
403    
404        // from here on we know that each sample consists of one value
405        FileWriter* fw = new FileWriter();
406        fw->openFile(filename, fileSize);
407        MPIBarrier();
408    
409        for (index_t y=0; y<myN1; y++) {
410            const int fileofs = (m_offset[0]+(m_offset[1]+y)*totalN0)*sizeof(ValueType);
411            ostringstream oss;
412    
413            for (index_t x=0; x<myN0; x++) {
414                const double* sample = _in->getSampleDataRO(y*myN0+x);
415                ValueType fvalue = static_cast<ValueType>(*sample);
416                if (byteOrder == BYTEORDER_NATIVE) {
417                    oss.write((char*)&fvalue, sizeof(fvalue));
418                } else {
419                    char* value = reinterpret_cast<char*>(&fvalue);
420                    oss.write(byte_swap32(value), sizeof(fvalue));
421                }
422            }
423            fw->writeAt(oss, fileofs);
424        }
425        fw->close();
426    }
427    
428  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
429  {  {
430  #if USE_SILO  #if USE_SILO
# Line 168  void Rectangle::dump(const string& fileN Line 487  void Rectangle::dump(const string& fileN
487      }      }
488      */      */
489    
490      boost::scoped_ptr<double> x(new double[m_N0]);      boost::scoped_ptr<double> x(new double[m_NN[0]]);
491      boost::scoped_ptr<double> y(new double[m_N1]);      boost::scoped_ptr<double> y(new double[m_NN[1]]);
492      double* coords[2] = { x.get(), y.get() };      double* coords[2] = { x.get(), y.get() };
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
493  #pragma omp parallel  #pragma omp parallel
494      {      {
495  #pragma omp for nowait  #pragma omp for nowait
496          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
497              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=getLocalCoordinate(i0, 0);
498          }          }
499  #pragma omp for nowait  #pragma omp for nowait
500          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
501              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=getLocalCoordinate(i1, 1);
502          }          }
503      }      }
504      IndexVector dims = getNumNodesPerDim();      int* dims = const_cast<int*>(getNumNodesPerDim());
505    
506      // write mesh      // write mesh
507      DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,
508              DB_COLLINEAR, NULL);              DB_COLLINEAR, NULL);
509    
510      // write node ids      // write node ids
511      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 2,      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 2,
512              NULL, 0, DB_INT, DB_NODECENT, NULL);              NULL, 0, DB_INT, DB_NODECENT, NULL);
513    
514      // write element ids      // write element ids
515      dims = getNumElementsPerDim();      dims = const_cast<int*>(getNumElementsPerDim());
516      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
517              &dims[0], 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);              dims, 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);
518    
519      // rank 0 writes multimesh and multivar      // rank 0 writes multimesh and multivar
520      if (m_mpiInfo->rank == 0) {      if (m_mpiInfo->rank == 0) {
# Line 289  bool Rectangle::ownSample(int fsType, in Line 606  bool Rectangle::ownSample(int fsType, in
606          case Elements:          case Elements:
607          case ReducedElements:          case ReducedElements:
608              // check ownership of element's bottom left node              // check ownership of element's bottom left node
609              return (m_dofMap[id%m_NE0+m_N0*(id/m_NE0)] < getNumDOF());              return (m_dofMap[id%m_NE[0]+m_NN[0]*(id/m_NE[0])] < getNumDOF());
610          case FaceElements:          case FaceElements:
611          case ReducedFaceElements:          case ReducedFaceElements:
612              {              {
613                  // determine which face the sample belongs to before                  // determine which face the sample belongs to before
614                  // checking ownership of corresponding element's first node                  // checking ownership of corresponding element's first node
                 const IndexVector faces = getNumFacesPerBoundary();  
615                  dim_t n=0;                  dim_t n=0;
616                  for (size_t i=0; i<faces.size(); i++) {                  for (size_t i=0; i<4; i++) {
617                      n+=faces[i];                      n+=m_faceCount[i];
618                      if (id<n) {                      if (id<n) {
619                          index_t k;                          index_t k;
620                          if (i==1)                          if (i==1)
621                              k=m_N0-2;                              k=m_NN[0]-2;
622                          else if (i==3)                          else if (i==3)
623                              k=m_N0*(m_N1-2);                              k=m_NN[0]*(m_NN[1]-2);
624                          else                          else
625                              k=0;                              k=0;
626                          // determine whether to move right or up                          // determine whether to move right or up
627                          const index_t delta=(i/2==0 ? m_N0 : 1);                          const index_t delta=(i/2==0 ? m_NN[0] : 1);
628                          return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());                          return (m_dofMap[k+(id-n+m_faceCount[i])*delta] < getNumDOF());
629                      }                      }
630                  }                  }
631                  return false;                  return false;
# Line 331  void Rectangle::setToNormal(escript::Dat Line 647  void Rectangle::setToNormal(escript::Dat
647          {          {
648              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
649  #pragma omp for nowait  #pragma omp for nowait
650                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
651                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
652                      // set vector at two quadrature points                      // set vector at two quadrature points
653                      *o++ = -1.;                      *o++ = -1.;
# Line 343  void Rectangle::setToNormal(escript::Dat Line 659  void Rectangle::setToNormal(escript::Dat
659    
660              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
661  #pragma omp for nowait  #pragma omp for nowait
662                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
663                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
664                      // set vector at two quadrature points                      // set vector at two quadrature points
665                      *o++ = 1.;                      *o++ = 1.;
# Line 355  void Rectangle::setToNormal(escript::Dat Line 671  void Rectangle::setToNormal(escript::Dat
671    
672              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
673  #pragma omp for nowait  #pragma omp for nowait
674                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
675                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
676                      // set vector at two quadrature points                      // set vector at two quadrature points
677                      *o++ = 0.;                      *o++ = 0.;
# Line 367  void Rectangle::setToNormal(escript::Dat Line 683  void Rectangle::setToNormal(escript::Dat
683    
684              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
685  #pragma omp for nowait  #pragma omp for nowait
686                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
687                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
688                      // set vector at two quadrature points                      // set vector at two quadrature points
689                      *o++ = 0.;                      *o++ = 0.;
# Line 383  void Rectangle::setToNormal(escript::Dat Line 699  void Rectangle::setToNormal(escript::Dat
699          {          {
700              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
701  #pragma omp for nowait  #pragma omp for nowait
702                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
703                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
704                      *o++ = -1.;                      *o++ = -1.;
705                      *o = 0.;                      *o = 0.;
# Line 392  void Rectangle::setToNormal(escript::Dat Line 708  void Rectangle::setToNormal(escript::Dat
708    
709              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
710  #pragma omp for nowait  #pragma omp for nowait
711                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
712                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
713                      *o++ = 1.;                      *o++ = 1.;
714                      *o = 0.;                      *o = 0.;
# Line 401  void Rectangle::setToNormal(escript::Dat Line 717  void Rectangle::setToNormal(escript::Dat
717    
718              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
719  #pragma omp for nowait  #pragma omp for nowait
720                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
721                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
722                      *o++ = 0.;                      *o++ = 0.;
723                      *o = -1.;                      *o = -1.;
# Line 410  void Rectangle::setToNormal(escript::Dat Line 726  void Rectangle::setToNormal(escript::Dat
726    
727              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
728  #pragma omp for nowait  #pragma omp for nowait
729                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
730                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
731                      *o++ = 0.;                      *o++ = 0.;
732                      *o = 1.;                      *o = 1.;
# Line 432  void Rectangle::setToSize(escript::Data& Line 748  void Rectangle::setToSize(escript::Data&
748              || out.getFunctionSpace().getTypeCode() == ReducedElements) {              || out.getFunctionSpace().getTypeCode() == ReducedElements) {
749          out.requireWrite();          out.requireWrite();
750          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
751          const double hSize=getFirstCoordAndSpacing(0).second;          const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]);
         const double vSize=getFirstCoordAndSpacing(1).second;  
         const double size=min(hSize,vSize);  
752  #pragma omp parallel for  #pragma omp parallel for
753          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
754              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 444  void Rectangle::setToSize(escript::Data& Line 758  void Rectangle::setToSize(escript::Data&
758              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
759          out.requireWrite();          out.requireWrite();
760          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
         const double hSize=getFirstCoordAndSpacing(0).second;  
         const double vSize=getFirstCoordAndSpacing(1).second;  
761  #pragma omp parallel  #pragma omp parallel
762          {          {
763              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
764  #pragma omp for nowait  #pragma omp for nowait
765                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
766                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
767                      fill(o, o+numQuad, vSize);                      fill(o, o+numQuad, m_dx[1]);
768                  }                  }
769              }              }
770    
771              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
772  #pragma omp for nowait  #pragma omp for nowait
773                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
774                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
775                      fill(o, o+numQuad, vSize);                      fill(o, o+numQuad, m_dx[1]);
776                  }                  }
777              }              }
778    
779              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
780  #pragma omp for nowait  #pragma omp for nowait
781                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
782                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
783                      fill(o, o+numQuad, hSize);                      fill(o, o+numQuad, m_dx[0]);
784                  }                  }
785              }              }
786    
787              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
788  #pragma omp for nowait  #pragma omp for nowait
789                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
790                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
791                      fill(o, o+numQuad, hSize);                      fill(o, o+numQuad, m_dx[0]);
792                  }                  }
793              }              }
794          } // end of parallel section          } // end of parallel section
# Line 489  void Rectangle::setToSize(escript::Data& Line 801  void Rectangle::setToSize(escript::Data&
801      }      }
802  }  }
803    
 Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  
                                                 bool reducedColOrder) const  
 {  
     /* FIXME: reduced  
     if (reducedRowOrder || reducedColOrder)  
         throw RipleyException("getPattern() not implemented for reduced order");  
     */  
     return m_pattern;  
 }  
   
804  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
805  {  {
806      RipleyDomain::Print_Mesh_Info(full);      RipleyDomain::Print_Mesh_Info(full);
# Line 506  void Rectangle::Print_Mesh_Info(const bo Line 808  void Rectangle::Print_Mesh_Info(const bo
808          cout << "     Id  Coordinates" << endl;          cout << "     Id  Coordinates" << endl;
809          cout.precision(15);          cout.precision(15);
810          cout.setf(ios::scientific, ios::floatfield);          cout.setf(ios::scientific, ios::floatfield);
         pair<double,double> xdx = getFirstCoordAndSpacing(0);  
         pair<double,double> ydy = getFirstCoordAndSpacing(1);  
811          for (index_t i=0; i < getNumNodes(); i++) {          for (index_t i=0; i < getNumNodes(); i++) {
812              cout << "  " << setw(5) << m_nodeId[i]              cout << "  " << setw(5) << m_nodeId[i]
813                  << "  " << xdx.first+(i%m_N0)*xdx.second                  << "  " << getLocalCoordinate(i%m_NN[0], 0)
814                  << "  " << ydy.first+(i/m_N0)*ydy.second << endl;                  << "  " << getLocalCoordinate(i/m_NN[0], 1) << endl;
815          }          }
816      }      }
817  }  }
818    
 IndexVector Rectangle::getNumNodesPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_N0);  
     ret.push_back(m_N1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumElementsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NE0);  
     ret.push_back(m_NE1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumFacesPerBoundary() const  
 {  
     IndexVector ret(4, 0);  
     //left  
     if (m_offset0==0)  
         ret[0]=m_NE1;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         ret[1]=m_NE1;  
     //bottom  
     if (m_offset1==0)  
         ret[2]=m_NE0;  
     //top  
     if (m_mpiInfo->rank/m_NX==m_NY-1)  
         ret[3]=m_NE0;  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumSubdivisionsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NX);  
     ret.push_back(m_NY);  
     return ret;  
 }  
   
 pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  
 {  
     if (dim==0) {  
         return pair<double,double>(m_x0+(m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);  
     } else if (dim==1) {  
         return pair<double,double>(m_y0+(m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);  
     }  
     throw RipleyException("getFirstCoordAndSpacing: invalid argument");  
 }  
   
 //protected  
 dim_t Rectangle::getNumDOF() const  
 {  
     return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;  
 }  
   
 //protected  
 dim_t Rectangle::getNumFaceElements() const  
 {  
     const IndexVector faces = getNumFacesPerBoundary();  
     dim_t n=0;  
     for (size_t i=0; i<faces.size(); i++)  
         n+=faces[i];  
     return n;  
 }  
819    
820  //protected  //protected
821  void Rectangle::assembleCoordinates(escript::Data& arg) const  void Rectangle::assembleCoordinates(escript::Data& arg) const
# Line 594  void Rectangle::assembleCoordinates(escr Line 827  void Rectangle::assembleCoordinates(escr
827      if (!numSamplesEqual(&x, 1, getNumNodes()))      if (!numSamplesEqual(&x, 1, getNumNodes()))
828          throw RipleyException("setToX: Illegal number of samples in Data object");          throw RipleyException("setToX: Illegal number of samples in Data object");
829    
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
830      arg.requireWrite();      arg.requireWrite();
831  #pragma omp parallel for  #pragma omp parallel for
832      for (dim_t i1 = 0; i1 < m_N1; i1++) {      for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
833          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
834              double* point = arg.getSampleDataRW(i0+m_N0*i1);              double* point = arg.getSampleDataRW(i0+m_NN[0]*i1);
835              point[0] = xdx.first+i0*xdx.second;              point[0] = getLocalCoordinate(i0, 0);
836              point[1] = ydy.first+i1*ydy.second;              point[1] = getLocalCoordinate(i1, 1);
837          }          }
838      }      }
839  }  }
# Line 611  void Rectangle::assembleCoordinates(escr Line 842  void Rectangle::assembleCoordinates(escr
842  void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const  void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
843  {  {
844      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
845      const double h0 = m_l0/m_gNE0;      const double cx0 = -1./m_dx[0];
846      const double h1 = m_l1/m_gNE1;      const double cx1 = -.78867513459481288225/m_dx[0];
847      const double cx0 = -1./h0;      const double cx2 = -.5/m_dx[0];
848      const double cx1 = -.78867513459481288225/h0;      const double cx3 = -.21132486540518711775/m_dx[0];
849      const double cx2 = -.5/h0;      const double cx4 = .21132486540518711775/m_dx[0];
850      const double cx3 = -.21132486540518711775/h0;      const double cx5 = .5/m_dx[0];
851      const double cx4 = .21132486540518711775/h0;      const double cx6 = .78867513459481288225/m_dx[0];
852      const double cx5 = .5/h0;      const double cx7 = 1./m_dx[0];
853      const double cx6 = .78867513459481288225/h0;      const double cy0 = -1./m_dx[1];
854      const double cx7 = 1./h0;      const double cy1 = -.78867513459481288225/m_dx[1];
855      const double cy0 = -1./h1;      const double cy2 = -.5/m_dx[1];
856      const double cy1 = -.78867513459481288225/h1;      const double cy3 = -.21132486540518711775/m_dx[1];
857      const double cy2 = -.5/h1;      const double cy4 = .21132486540518711775/m_dx[1];
858      const double cy3 = -.21132486540518711775/h1;      const double cy5 = .5/m_dx[1];
859      const double cy4 = .21132486540518711775/h1;      const double cy6 = .78867513459481288225/m_dx[1];
860      const double cy5 = .5/h1;      const double cy7 = 1./m_dx[1];
     const double cy6 = .78867513459481288225/h1;  
     const double cy7 = 1./h1;  
861    
862      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
863          out.requireWrite();          out.requireWrite();
864  #pragma omp parallel for  #pragma omp parallel
865          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
866              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
867                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
868                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
869                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
870                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
871                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
872                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
873                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
874                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
875                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
876                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
877                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
878                      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) {
879                      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;
880                      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;
881                  } /* 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;
882              } /* 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;
883          } /* 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;
884                            o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
885                            o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
886                            o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
887                        } // end of component loop i
888                    } // end of k0 loop
889                } // end of k1 loop
890            } // end of parallel section
891      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
892          out.requireWrite();          out.requireWrite();
893  #pragma omp parallel for  #pragma omp parallel
894          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
895              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
896                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
897                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
898                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
899                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
900                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
901                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
902                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
903                      o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
904                  } /* end of component loop i */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
905              } /* end of k0 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
906          } /* end of k1 loop */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
907                        for (index_t i=0; i < numComp; ++i) {
908                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
909                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
910                        } // end of component loop i
911                    } // end of k0 loop
912                } // end of k1 loop
913            } // end of parallel section
914      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
915          out.requireWrite();          out.requireWrite();
916  #pragma omp parallel  #pragma omp parallel
917          {          {
918                vector<double> f_00(numComp);
919                vector<double> f_01(numComp);
920                vector<double> f_10(numComp);
921                vector<double> f_11(numComp);
922              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
923  #pragma omp for nowait  #pragma omp for nowait
924                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
925                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
926                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
927                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
928                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
929                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
930                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
931                          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;
932                          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;
933                          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;
934                          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;
935                      } /* end of component loop i */                      } // end of component loop i
936                  } /* end of k1 loop */                  } // end of k1 loop
937              } /* end of face 0 */              } // end of face 0
938              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
939  #pragma omp for nowait  #pragma omp for nowait
940                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
941                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
942                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
943                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
944                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
945                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
946                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
947                          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;
948                          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;
949                          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;
950                          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;
951                      } /* end of component loop i */                      } // end of component loop i
952                  } /* end of k1 loop */                  } // end of k1 loop
953              } /* end of face 1 */              } // end of face 1
954              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
955  #pragma omp for nowait  #pragma omp for nowait
956                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
957                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
958                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
959                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
960                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
961                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
962                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
963                          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;
964                          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;
965                          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;
966                          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;
967                      } /* end of component loop i */                      } // end of component loop i
968                  } /* end of k0 loop */                  } // end of k0 loop
969              } /* end of face 2 */              } // end of face 2
970              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
971  #pragma omp for nowait  #pragma omp for nowait
972                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
973                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
974                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
975                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
976                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
977                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
978                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
979                          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;
980                          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;
981                          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;
982                          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;
983                      } /* end of component loop i */                      } // end of component loop i
984                  } /* end of k0 loop */                  } // end of k0 loop
985              } /* end of face 3 */              } // end of face 3
986          } // end of parallel section          } // end of parallel section
987    
988      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
989          out.requireWrite();          out.requireWrite();
990  #pragma omp parallel  #pragma omp parallel
991          {          {
992                vector<double> f_00(numComp);
993                vector<double> f_01(numComp);
994                vector<double> f_10(numComp);
995                vector<double> f_11(numComp);
996              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
997  #pragma omp for nowait  #pragma omp for nowait
998                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
999                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1000                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1001                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_NN[0])), numComp*sizeof(double));
1002                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_NN[0])), numComp*sizeof(double));
1003                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1004                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1005                          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]);
1006                          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;
1007                      } /* end of component loop i */                      } // end of component loop i
1008                  } /* end of k1 loop */                  } // end of k1 loop
1009              } /* end of face 0 */              } // end of face 0
1010              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1011  #pragma omp for nowait  #pragma omp for nowait
1012                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1013                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1, m_NN[0])), numComp*sizeof(double));
1014                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_NN[0]-2,k1+1, m_NN[0])), numComp*sizeof(double));
1015                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1016                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1017                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1018                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1019                          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]);
1020                          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;
1021                      } /* end of component loop i */                      } // end of component loop i
1022                  } /* end of k1 loop */                  } // end of k1 loop
1023              } /* end of face 1 */              } // end of face 1
1024              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1025  #pragma omp for nowait  #pragma omp for nowait
1026                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1027                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1028                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_NN[0])), numComp*sizeof(double));
1029                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1030                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_NN[0])), numComp*sizeof(double));
1031                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1032                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1033                          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;
1034                          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]);
1035                      } /* end of component loop i */                      } // end of component loop i
1036                  } /* end of k0 loop */                  } // end of k0 loop
1037              } /* end of face 2 */              } // end of face 2
1038              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1039  #pragma omp for nowait  #pragma omp for nowait
1040                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1041                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1042                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1043                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1044                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1045                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1046                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1047                          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;
1048                          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]);
1049                      } /* end of component loop i */                      } // end of component loop i
1050                  } /* end of k0 loop */                  } // end of k0 loop
1051              } /* end of face 3 */              } // end of face 3
1052          } // end of parallel section          } // end of parallel section
1053      }      }
1054  }  }
# Line 805  void Rectangle::assembleGradient(escript Line 1057  void Rectangle::assembleGradient(escript
1057  void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const  void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1058  {  {
1059      const dim_t numComp = arg.getDataPointSize();      const dim_t numComp = arg.getDataPointSize();
1060      const double h0 = m_l0/m_gNE0;      const index_t left = (m_offset[0]==0 ? 0 : 1);
1061      const double h1 = m_l1/m_gNE1;      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1062      const index_t left = (m_offset0==0 ? 0 : 1);      const int fs=arg.getFunctionSpace().getTypeCode();
1063      const index_t bottom = (m_offset1==0 ? 0 : 1);      if (fs == Elements && arg.actsExpanded()) {
     if (arg.getFunctionSpace().getTypeCode() == Elements) {  
         const double w = h0*h1/4.;  
1064  #pragma omp parallel  #pragma omp parallel
1065          {          {
1066              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1067                const double w = m_dx[0]*m_dx[1]/4.;
1068  #pragma omp for nowait  #pragma omp for nowait
1069              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1070                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1071                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1072                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1073                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1074                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1075                          const double f2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
1076                          const double f3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
1077                          int_local[i]+=(f0+f1+f2+f3)*w;                          int_local[i]+=(f0+f1+f2+f3)*w;
1078                      }  /* end of component loop i */                      }  // end of component loop i
1079                  } /* end of k0 loop */                  } // end of k0 loop
1080              } /* end of k1 loop */              } // end of k1 loop
   
1081  #pragma omp critical  #pragma omp critical
1082              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1083                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1084          } // end of parallel section          } // end of parallel section
1085      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
1086          const double w = h0*h1;      } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1087            const double w = m_dx[0]*m_dx[1];
1088  #pragma omp parallel  #pragma omp parallel
1089          {          {
1090              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1091  #pragma omp for nowait  #pragma omp for nowait
1092              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1093                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1094                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1095                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1096                          int_local[i]+=f[i]*w;                          int_local[i]+=f[i]*w;
1097                      }  /* end of component loop i */                      }
1098                  } /* end of k0 loop */                  }
1099              } /* end of k1 loop */              }
   
1100  #pragma omp critical  #pragma omp critical
1101              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1102                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1103          } // end of parallel section          } // end of parallel section
1104      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
1105          const double w0 = h0/2.;      } else if (fs == FaceElements && arg.actsExpanded()) {
         const double w1 = h1/2.;  
1106  #pragma omp parallel  #pragma omp parallel
1107          {          {
1108              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1109                const double w0 = m_dx[0]/2.;
1110                const double w1 = m_dx[1]/2.;
1111              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1112  #pragma omp for nowait  #pragma omp for nowait
1113                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1114                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1115                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1116                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1117                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1118                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1119                      }  /* end of component loop i */                      }  // end of component loop i
1120                  } /* end of k1 loop */                  } // end of k1 loop
1121              }              }
1122    
1123              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1124  #pragma omp for nowait  #pragma omp for nowait
1125                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1126                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1127                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1128                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1129                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1130                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1131                      }  /* end of component loop i */                      }  // end of component loop i
1132                  } /* end of k1 loop */                  } // end of k1 loop
1133              }              }
1134    
1135              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1136  #pragma omp for nowait  #pragma omp for nowait
1137                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1138                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1139                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1140                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1141                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1142                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1143                      }  /* end of component loop i */                      }  // end of component loop i
1144                  } /* end of k0 loop */                  } // end of k0 loop
1145              }              }
1146    
1147              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1148  #pragma omp for nowait  #pragma omp for nowait
1149                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1150                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1151                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1152                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1153                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1154                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1155                      }  /* end of component loop i */                      }  // end of component loop i
1156                  } /* end of k0 loop */                  } // end of k0 loop
1157              }              }
   
1158  #pragma omp critical  #pragma omp critical
1159              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1160                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1161          } // end of parallel section          } // end of parallel section
1162      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
1163        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1164  #pragma omp parallel  #pragma omp parallel
1165          {          {
1166              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1167              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1168  #pragma omp for nowait  #pragma omp for nowait
1169                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1170                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1171                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1172                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*m_dx[1];
1173                      }  /* end of component loop i */                      }
1174                  } /* end of k1 loop */                  }
1175              }              }
1176    
1177              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1178  #pragma omp for nowait  #pragma omp for nowait
1179                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1180                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1181                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1182                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*m_dx[1];
1183                      }  /* end of component loop i */                      }
1184                  } /* end of k1 loop */                  }
1185              }              }
1186    
1187              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1188  #pragma omp for nowait  #pragma omp for nowait
1189                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1190                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1191                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1192                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*m_dx[0];
1193                      }  /* end of component loop i */                      }
1194                  } /* end of k0 loop */                  }
1195              }              }
1196    
1197              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1198  #pragma omp for nowait  #pragma omp for nowait
1199                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1200                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1201                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1202                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*m_dx[0];
1203                      }  /* end of component loop i */                      }
1204                  } /* end of k0 loop */                  }
1205              }              }
1206    
1207  #pragma omp critical  #pragma omp critical
1208              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1209                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1210          } // end of parallel section          } // end of parallel section
1211      }      } // function space selector
1212  }  }
1213    
1214  //protected  //protected
1215  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
1216  {  {
1217      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1218      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1219      const int x=node%nDOF0;      const int x=node%nDOF0;
1220      const int y=node/nDOF0;      const int y=node/nDOF0;
1221      dim_t num=0;      dim_t num=0;
# Line 994  void Rectangle::nodesToDOF(escript::Data Line 1245  void Rectangle::nodesToDOF(escript::Data
1245      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1246      out.requireWrite();      out.requireWrite();
1247    
1248      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1249      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1250      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1251      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1252  #pragma omp parallel for  #pragma omp parallel for
1253      for (index_t i=0; i<nDOF1; i++) {      for (index_t i=0; i<nDOF1; i++) {
1254          for (index_t j=0; j<nDOF0; j++) {          for (index_t j=0; j<nDOF0; j++) {
1255              const index_t n=j+left+(i+bottom)*m_N0;              const index_t n=j+left+(i+bottom)*m_NN[0];
1256              const double* src=in.getSampleDataRO(n);              const double* src=in.getSampleDataRO(n);
1257              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
1258          }          }
# Line 1027  void Rectangle::dofToNodes(escript::Data Line 1278  void Rectangle::dofToNodes(escript::Data
1278                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1279          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1280      }      }
1281        Paso_Coupler_free(coupler);
1282  }  }
1283    
1284  //private  //private
1285  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
1286  {  {
1287      // identifiers are ordered from left to right, bottom to top globablly.      // degrees of freedom are numbered from left to right, bottom to top in
1288        // each rank, continuing on the next rank (ranks also go left-right,
1289        // bottom-top).
1290        // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
1291        // helps when writing out data rank after rank.
1292    
1293      // build node distribution vector first.      // build node distribution vector first.
1294      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
1295        // constant for all ranks in this implementation
1296      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1297      const dim_t numDOF=getNumDOF();      const dim_t numDOF=getNumDOF();
1298      for (dim_t k=1; k<m_mpiInfo->size; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
# Line 1045  void Rectangle::populateSampleIds() Line 1302  void Rectangle::populateSampleIds()
1302      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
1303      m_dofId.resize(numDOF);      m_dofId.resize(numDOF);
1304      m_elementId.resize(getNumElements());      m_elementId.resize(getNumElements());
1305    
1306        // populate face element counts
1307        //left
1308        if (m_offset[0]==0)
1309            m_faceCount[0]=m_NE[1];
1310        else
1311            m_faceCount[0]=0;
1312        //right
1313        if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
1314            m_faceCount[1]=m_NE[1];
1315        else
1316            m_faceCount[1]=0;
1317        //bottom
1318        if (m_offset[1]==0)
1319            m_faceCount[2]=m_NE[0];
1320        else
1321            m_faceCount[2]=0;
1322        //top
1323        if (m_mpiInfo->rank/m_NX[0]==m_NX[1]-1)
1324            m_faceCount[3]=m_NE[0];
1325        else
1326            m_faceCount[3]=0;
1327    
1328      m_faceId.resize(getNumFaceElements());      m_faceId.resize(getNumFaceElements());
1329    
1330        const index_t left = (m_offset[0]==0 ? 0 : 1);
1331        const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1332        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1333        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1334    
1335    #define globalNodeId(x,y) \
1336        ((m_offset[0]+x)/nDOF0)*nDOF0*nDOF1+(m_offset[0]+x)%nDOF0 \
1337        + ((m_offset[1]+y)/nDOF1)*nDOF0*nDOF1*m_NX[0]+((m_offset[1]+y)%nDOF1)*nDOF0
1338    
1339        // set corner id's outside the parallel region
1340        m_nodeId[0] = globalNodeId(0, 0);
1341        m_nodeId[m_NN[0]-1] = globalNodeId(m_NN[0]-1, 0);
1342        m_nodeId[m_NN[0]*(m_NN[1]-1)] = globalNodeId(0, m_NN[1]-1);
1343        m_nodeId[m_NN[0]*m_NN[1]-1] = globalNodeId(m_NN[0]-1,m_NN[1]-1);
1344    #undef globalNodeId
1345    
1346  #pragma omp parallel  #pragma omp parallel
1347      {      {
1348          // nodes          // populate degrees of freedom and own nodes (identical id)
1349  #pragma omp for nowait  #pragma omp for nowait
1350          for (dim_t i1=0; i1<m_N1; i1++) {          for (dim_t i=0; i<nDOF1; i++) {
1351              for (dim_t i0=0; i0<m_N0; i0++) {              for (dim_t j=0; j<nDOF0; j++) {
1352                  m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;                  const index_t nodeIdx=j+left+(i+bottom)*m_NN[0];
1353                    const index_t dofIdx=j+i*nDOF0;
1354                    m_dofId[dofIdx] = m_nodeId[nodeIdx]
1355                        = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
1356              }              }
1357          }          }
1358    
1359          // degrees of freedom          // populate the rest of the nodes (shared with other ranks)
1360            if (m_faceCount[0]==0) { // left column
1361  #pragma omp for nowait  #pragma omp for nowait
1362          for (dim_t k=0; k<numDOF; k++)              for (dim_t i=0; i<nDOF1; i++) {
1363              m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;                  const index_t nodeIdx=(i+bottom)*m_NN[0];
1364                    const index_t dofId=(i+1)*nDOF0-1;
1365                    m_nodeId[nodeIdx]
1366                        = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
1367                }
1368            }
1369            if (m_faceCount[1]==0) { // right column
1370    #pragma omp for nowait
1371                for (dim_t i=0; i<nDOF1; i++) {
1372                    const index_t nodeIdx=(i+bottom+1)*m_NN[0]-1;
1373                    const index_t dofId=i*nDOF0;
1374                    m_nodeId[nodeIdx]
1375                        = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
1376                }
1377            }
1378            if (m_faceCount[2]==0) { // bottom row
1379    #pragma omp for nowait
1380                for (dim_t i=0; i<nDOF0; i++) {
1381                    const index_t nodeIdx=i+left;
1382                    const index_t dofId=nDOF0*(nDOF1-1)+i;
1383                    m_nodeId[nodeIdx]
1384                        = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
1385                }
1386            }
1387            if (m_faceCount[3]==0) { // top row
1388    #pragma omp for nowait
1389                for (dim_t i=0; i<nDOF0; i++) {
1390                    const index_t nodeIdx=m_NN[0]*(m_NN[1]-1)+i+left;
1391                    const index_t dofId=i;
1392                    m_nodeId[nodeIdx]
1393                        = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
1394                }
1395            }
1396    
1397          // elements          // populate element id's
1398  #pragma omp for nowait  #pragma omp for nowait
1399          for (dim_t i1=0; i1<m_NE1; i1++) {          for (dim_t i1=0; i1<m_NE[1]; i1++) {
1400              for (dim_t i0=0; i0<m_NE0; i0++) {              for (dim_t i0=0; i0<m_NE[0]; i0++) {
1401                  m_elementId[i0+i1*m_NE0]=(m_offset1+i1)*m_gNE0+m_offset0+i0;                  m_elementId[i0+i1*m_NE[0]]=(m_offset[1]+i1)*m_gNE[0]+m_offset[0]+i0;
1402              }              }
1403          }          }
1404    
# Line 1083  void Rectangle::populateSampleIds() Line 1415  void Rectangle::populateSampleIds()
1415      updateTagsInUse(Elements);      updateTagsInUse(Elements);
1416    
1417      // generate face offset vector and set face tags      // generate face offset vector and set face tags
     const IndexVector facesPerEdge = getNumFacesPerBoundary();  
1418      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
1419      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
1420      m_faceOffset.assign(facesPerEdge.size(), -1);      m_faceOffset.assign(4, -1);
1421      m_faceTags.clear();      m_faceTags.clear();
1422      index_t offset=0;      index_t offset=0;
1423      for (size_t i=0; i<facesPerEdge.size(); i++) {      for (size_t i=0; i<4; i++) {
1424          if (facesPerEdge[i]>0) {          if (m_faceCount[i]>0) {
1425              m_faceOffset[i]=offset;              m_faceOffset[i]=offset;
1426              offset+=facesPerEdge[i];              offset+=m_faceCount[i];
1427              m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);              m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
1428          }          }
1429      }      }
1430      setTagMap("left", LEFT);      setTagMap("left", LEFT);
# Line 1106  void Rectangle::populateSampleIds() Line 1437  void Rectangle::populateSampleIds()
1437  //private  //private
1438  void Rectangle::createPattern()  void Rectangle::createPattern()
1439  {  {
1440      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1441      const dim_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1442      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1443      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1444    
1445      // populate node->DOF mapping with own degrees of freedom.      // populate node->DOF mapping with own degrees of freedom.
1446      // The rest is assigned in the loop further down      // The rest is assigned in the loop further down
# Line 1117  void Rectangle::createPattern() Line 1448  void Rectangle::createPattern()
1448  #pragma omp parallel for  #pragma omp parallel for
1449      for (index_t i=bottom; i<bottom+nDOF1; i++) {      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1450          for (index_t j=left; j<left+nDOF0; j++) {          for (index_t j=left; j<left+nDOF0; j++) {
1451              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;              m_dofMap[i*m_NN[0]+j]=(i-bottom)*nDOF0+j-left;
1452          }          }
1453      }      }
1454    
# Line 1130  void Rectangle::createPattern() Line 1461  void Rectangle::createPattern()
1461      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
1462      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
1463      int numShared=0;      int numShared=0;
1464      const int x=m_mpiInfo->rank%m_NX;      const int x=m_mpiInfo->rank%m_NX[0];
1465      const int y=m_mpiInfo->rank/m_NX;      const int y=m_mpiInfo->rank/m_NX[0];
1466      for (int i1=-1; i1<2; i1++) {      for (int i1=-1; i1<2; i1++) {
1467          for (int i0=-1; i0<2; i0++) {          for (int i0=-1; i0<2; i0++) {
1468              // skip this rank              // skip this rank
# Line 1140  void Rectangle::createPattern() Line 1471  void Rectangle::createPattern()
1471              // location of neighbour rank              // location of neighbour rank
1472              const int nx=x+i0;              const int nx=x+i0;
1473              const int ny=y+i1;              const int ny=y+i1;
1474              if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {              if (nx>=0 && ny>=0 && nx<m_NX[0] && ny<m_NX[1]) {
1475                  neighbour.push_back(ny*m_NX+nx);                  neighbour.push_back(ny*m_NX[0]+nx);
1476                  if (i0==0) {                  if (i0==0) {
1477                      // sharing top or bottom edge                      // sharing top or bottom edge
1478                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1479                      const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);                      const int firstNode=(i1==-1 ? left : m_NN[0]*(m_NN[1]-1)+left);
1480                      offsetInShared.push_back(offsetInShared.back()+nDOF0);                      offsetInShared.push_back(offsetInShared.back()+nDOF0);
1481                      for (dim_t i=0; i<nDOF0; i++, numShared++) {                      for (dim_t i=0; i<nDOF0; i++, numShared++) {
1482                          sendShared.push_back(firstDOF+i);                          sendShared.push_back(firstDOF+i);
# Line 1160  void Rectangle::createPattern() Line 1491  void Rectangle::createPattern()
1491                  } else if (i1==0) {                  } else if (i1==0) {
1492                      // sharing left or right edge                      // sharing left or right edge
1493                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1494                      const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);                      const int firstNode=(i0==-1 ? bottom*m_NN[0] : (bottom+1)*m_NN[0]-1);
1495                      offsetInShared.push_back(offsetInShared.back()+nDOF1);                      offsetInShared.push_back(offsetInShared.back()+nDOF1);
1496                      for (dim_t i=0; i<nDOF1; i++, numShared++) {                      for (dim_t i=0; i<nDOF1; i++, numShared++) {
1497                          sendShared.push_back(firstDOF+i*nDOF0);                          sendShared.push_back(firstDOF+i*nDOF0);
# Line 1170  void Rectangle::createPattern() Line 1501  void Rectangle::createPattern()
1501                          colIndices[firstDOF+i*nDOF0].push_back(numShared);                          colIndices[firstDOF+i*nDOF0].push_back(numShared);
1502                          if (i<nDOF1-1)                          if (i<nDOF1-1)
1503                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1504                          m_dofMap[firstNode+i*m_N0]=numDOF+numShared;                          m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
1505                      }                      }
1506                  } else {                  } else {
1507                      // sharing a node                      // sharing a node
1508                      const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);                      const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1509                      const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);                      const int node=(i0+1)/2*(m_NN[0]-1)+(i1+1)/2*m_NN[0]*(m_NN[1]-1);
1510                      offsetInShared.push_back(offsetInShared.back()+1);                      offsetInShared.push_back(offsetInShared.back()+1);
1511                      sendShared.push_back(dof);                      sendShared.push_back(dof);
1512                      recvShared.push_back(numDOF+numShared);                      recvShared.push_back(numDOF+numShared);
# Line 1285  void Rectangle::addToMatrixAndRHS(Paso_S Line 1616  void Rectangle::addToMatrixAndRHS(Paso_S
1616      IndexVector rowIndex;      IndexVector rowIndex;
1617      rowIndex.push_back(m_dofMap[firstNode]);      rowIndex.push_back(m_dofMap[firstNode]);
1618      rowIndex.push_back(m_dofMap[firstNode+1]);      rowIndex.push_back(m_dofMap[firstNode+1]);
1619      rowIndex.push_back(m_dofMap[firstNode+m_N0]);      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]]);
1620      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]+1]);
1621      if (addF) {      if (addF) {
1622          double *F_p=F.getSampleDataRW(0);          double *F_p=F.getSampleDataRW(0);
1623          for (index_t i=0; i<rowIndex.size(); i++) {          for (index_t i=0; i<rowIndex.size(); i++) {
# Line 1309  void Rectangle::interpolateNodesOnElemen Line 1640  void Rectangle::interpolateNodesOnElemen
1640      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1641      if (reduced) {      if (reduced) {
1642          out.requireWrite();          out.requireWrite();
1643          const double c0 = .25;          const double c0 = 0.25;
1644  #pragma omp parallel for  #pragma omp parallel
1645          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1646              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1647                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_01(numComp);
1648                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1649                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              vector<double> f_11(numComp);
1650                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1651                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1652                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1653                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
1654                  } /* end of component loop i */                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
1655              } /* end of k0 loop */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
1656          } /* end of k1 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
1657                        double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1658                        for (index_t i=0; i < numComp; ++i) {
1659                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1660                        } /* end of component loop i */
1661                    } /* end of k0 loop */
1662                } /* end of k1 loop */
1663            } /* end of parallel section */
1664      } else {      } else {
1665          out.requireWrite();          out.requireWrite();
1666          const double c0 = .16666666666666666667;          const double c0 = 0.16666666666666666667;
1667          const double c1 = .044658198738520451079;          const double c1 = 0.044658198738520451079;
1668          const double c2 = .62200846792814621559;          const double c2 = 0.62200846792814621559;
1669  #pragma omp parallel for  #pragma omp parallel
1670          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1671              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1672                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
1673                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1674                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
1675                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1676                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1677                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1678                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_NN[0])), numComp*sizeof(double));
1679                      o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
1680                      o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
1681                      o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
1682                  } /* end of component loop i */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1683              } /* end of k0 loop */                      for (index_t i=0; i < numComp; ++i) {
1684          } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1685                            o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1686                            o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1687                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1688                        } /* end of component loop i */
1689                    } /* end of k0 loop */
1690                } /* end of k1 loop */
1691            } /* end of parallel section */
1692      }      }
1693  }  }
1694    
# Line 1354  void Rectangle::interpolateNodesOnFaces( Line 1699  void Rectangle::interpolateNodesOnFaces(
1699      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1700      if (reduced) {      if (reduced) {
1701          out.requireWrite();          out.requireWrite();
1702          const double c0 = .5;          const double c0 = 0.5;
1703  #pragma omp parallel  #pragma omp parallel
1704          {          {
1705                vector<double> f_00(numComp);
1706                vector<double> f_01(numComp);
1707                vector<double> f_10(numComp);
1708                vector<double> f_11(numComp);
1709              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1710  #pragma omp for nowait  #pragma omp for nowait
1711                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1712                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1713                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1714                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1715                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1716                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
# Line 1370  void Rectangle::interpolateNodesOnFaces( Line 1719  void Rectangle::interpolateNodesOnFaces(
1719              } /* end of face 0 */              } /* end of face 0 */
1720              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1721  #pragma omp for nowait  #pragma omp for nowait
1722                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1723                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1724                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1725                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1726                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1727                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
# Line 1381  void Rectangle::interpolateNodesOnFaces( Line 1730  void Rectangle::interpolateNodesOnFaces(
1730              } /* end of face 1 */              } /* end of face 1 */
1731              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1732  #pragma omp for nowait  #pragma omp for nowait
1733                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1734                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1735                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1736                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1737                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1738                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
# Line 1392  void Rectangle::interpolateNodesOnFaces( Line 1741  void Rectangle::interpolateNodesOnFaces(
1741              } /* end of face 2 */              } /* end of face 2 */
1742              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1743  #pragma omp for nowait  #pragma omp for nowait
1744                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1745                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1746                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1747                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1748                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1749                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1750                      } /* end of component loop i */                      } /* end of component loop i */
1751                  } /* end of k0 loop */                  } /* end of k0 loop */
1752              } /* end of face 3 */              } /* end of face 3 */
1753          } // end of parallel section          } /* end of parallel section */
1754      } else {      } else {
1755          out.requireWrite();          out.requireWrite();
1756          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1757          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1758  #pragma omp parallel  #pragma omp parallel
1759          {          {
1760                vector<double> f_00(numComp);
1761                vector<double> f_01(numComp);
1762                vector<double> f_10(numComp);
1763                vector<double> f_11(numComp);
1764              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1765  #pragma omp for nowait  #pragma omp for nowait
1766                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1767                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_NN[0])), numComp*sizeof(double));
1768                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_NN[0])), numComp*sizeof(double));
1769                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1770                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1771                          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];
1772                          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];
1773                      } /* end of component loop i */                      } /* end of component loop i */
1774                  } /* end of k1 loop */                  } /* end of k1 loop */
1775              } /* end of face 0 */              } /* end of face 0 */
1776              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1777  #pragma omp for nowait  #pragma omp for nowait
1778                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1779                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1, m_NN[0])), numComp*sizeof(double));
1780                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_NN[0]-1,k1+1, m_NN[0])), numComp*sizeof(double));
1781                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1782                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1783                          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];
1784                          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];
1785                      } /* end of component loop i */                      } /* end of component loop i */
1786                  } /* end of k1 loop */                  } /* end of k1 loop */
1787              } /* end of face 1 */              } /* end of face 1 */
1788              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1789  #pragma omp for nowait  #pragma omp for nowait
1790                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1791                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_NN[0])), numComp*sizeof(double));
1792                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_NN[0])), numComp*sizeof(double));
1793                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1794                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1795                          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];
1796                          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];
1797                      } /* end of component loop i */                      } /* end of component loop i */
1798                  } /* end of k0 loop */                  } /* end of k0 loop */
1799              } /* end of face 2 */              } /* end of face 2 */
1800              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1801  #pragma omp for nowait  #pragma omp for nowait
1802                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1803                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1804                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1805                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1806                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1807                          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];
1808                          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];
1809                      } /* end of component loop i */                      } /* end of component loop i */
1810                  } /* end of k0 loop */                  } /* end of k0 loop */
1811              } /* end of face 3 */              } /* end of face 3 */
1812          } // end of parallel section          } /* end of parallel section */
1813      }      }
1814  }  }
1815    
# Line 1466  void Rectangle::assemblePDESingle(Paso_S Line 1819  void Rectangle::assemblePDESingle(Paso_S
1819          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
1820          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
1821  {  {
1822      const double h0 = m_l0/m_gNE0;      const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];
     const double h1 = m_l1/m_gNE1;  
     const double w0 = -0.1555021169820365539*h1/h0;  
1823      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
1824      const double w2 = -0.15550211698203655390;      const double w2 = -0.15550211698203655390;
1825      const double w3 = 0.041666666666666666667*h0/h1;      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];
1826      const double w4 = 0.15550211698203655390;      const double w4 = 0.15550211698203655390;
1827      const double w5 = -0.041666666666666666667;      const double w5 = -0.041666666666666666667;
1828      const double w6 = -0.01116454968463011277*h1/h0;      const double w6 = -0.01116454968463011277*m_dx[1]/m_dx[0];
1829      const double w7 = 0.011164549684630112770;      const double w7 = 0.011164549684630112770;
1830      const double w8 = -0.011164549684630112770;      const double w8 = -0.011164549684630112770;
1831      const double w9 = -0.041666666666666666667*h1/h0;      const double w9 = -0.041666666666666666667*m_dx[1]/m_dx[0];
1832      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.041666666666666666667*m_dx[0]/m_dx[1];
1833      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = 0.1555021169820365539*m_dx[1]/m_dx[0];
1834      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = 0.1555021169820365539*m_dx[0]/m_dx[1];
1835      const double w13 = 0.01116454968463011277*h0/h1;      const double w13 = 0.01116454968463011277*m_dx[0]/m_dx[1];
1836      const double w14 = 0.01116454968463011277*h1/h0;      const double w14 = 0.01116454968463011277*m_dx[1]/m_dx[0];
1837      const double w15 = 0.041666666666666666667*h1/h0;      const double w15 = 0.041666666666666666667*m_dx[1]/m_dx[0];
1838      const double w16 = -0.01116454968463011277*h0/h1;      const double w16 = -0.01116454968463011277*m_dx[0]/m_dx[1];
1839      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*m_dx[0]/m_dx[1];
1840      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*m_dx[1]/m_dx[0];
1841      const double w19 = 0.25;      const double w19 = 0.25;
1842      const double w20 = -0.25;      const double w20 = -0.25;
1843      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*m_dx[0]/m_dx[1];
1844      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*m_dx[1]/m_dx[0];
1845      const double w23 = -0.16666666666666666667*h0/h1;      const double w23 = -0.16666666666666666667*m_dx[0]/m_dx[1];
1846      const double w24 = 0.33333333333333333333*h1/h0;      const double w24 = 0.33333333333333333333*m_dx[1]/m_dx[0];
1847      const double w25 = 0.33333333333333333333*h0/h1;      const double w25 = 0.33333333333333333333*m_dx[0]/m_dx[1];
1848      const double w26 = 0.16666666666666666667*h1/h0;      const double w26 = 0.16666666666666666667*m_dx[1]/m_dx[0];
1849      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = -0.33333333333333333333*m_dx[0]/m_dx[1];
1850      const double w28 = -0.032861463941450536761*h1;      const double w28 = -0.032861463941450536761*m_dx[1];
1851      const double w29 = -0.032861463941450536761*h0;      const double w29 = -0.032861463941450536761*m_dx[0];
1852      const double w30 = -0.12264065304058601714*h1;      const double w30 = -0.12264065304058601714*m_dx[1];
1853      const double w31 = -0.0023593469594139828636*h1;      const double w31 = -0.0023593469594139828636*m_dx[1];
1854      const double w32 = -0.008805202725216129906*h0;      const double w32 = -0.008805202725216129906*m_dx[0];
1855      const double w33 = -0.008805202725216129906*h1;      const double w33 = -0.008805202725216129906*m_dx[1];
1856      const double w34 = 0.032861463941450536761*h1;      const double w34 = 0.032861463941450536761*m_dx[1];
1857      const double w35 = 0.008805202725216129906*h1;      const double w35 = 0.008805202725216129906*m_dx[1];
1858      const double w36 = 0.008805202725216129906*h0;      const double w36 = 0.008805202725216129906*m_dx[0];
1859      const double w37 = 0.0023593469594139828636*h1;      const double w37 = 0.0023593469594139828636*m_dx[1];
1860      const double w38 = 0.12264065304058601714*h1;      const double w38 = 0.12264065304058601714*m_dx[1];
1861      const double w39 = 0.032861463941450536761*h0;      const double w39 = 0.032861463941450536761*m_dx[0];
1862      const double w40 = -0.12264065304058601714*h0;      const double w40 = -0.12264065304058601714*m_dx[0];
1863      const double w41 = -0.0023593469594139828636*h0;      const double w41 = -0.0023593469594139828636*m_dx[0];
1864      const double w42 = 0.0023593469594139828636*h0;      const double w42 = 0.0023593469594139828636*m_dx[0];
1865      const double w43 = 0.12264065304058601714*h0;      const double w43 = 0.12264065304058601714*m_dx[0];
1866      const double w44 = -0.16666666666666666667*h1;      const double w44 = -0.16666666666666666667*m_dx[1];
1867      const double w45 = -0.083333333333333333333*h0;      const double w45 = -0.083333333333333333333*m_dx[0];
1868      const double w46 = 0.083333333333333333333*h1;      const double w46 = 0.083333333333333333333*m_dx[1];
1869      const double w47 = 0.16666666666666666667*h1;      const double w47 = 0.16666666666666666667*m_dx[1];
1870      const double w48 = 0.083333333333333333333*h0;      const double w48 = 0.083333333333333333333*m_dx[0];
1871      const double w49 = -0.16666666666666666667*h0;      const double w49 = -0.16666666666666666667*m_dx[0];
1872      const double w50 = 0.16666666666666666667*h0;      const double w50 = 0.16666666666666666667*m_dx[0];
1873      const double w51 = -0.083333333333333333333*h1;      const double w51 = -0.083333333333333333333*m_dx[1];
1874      const double w52 = 0.025917019497006092316*h0*h1;      const double w52 = 0.025917019497006092316*m_dx[0]*m_dx[1];
1875      const double w53 = 0.0018607582807716854616*h0*h1;      const double w53 = 0.0018607582807716854616*m_dx[0]*m_dx[1];
1876      const double w54 = 0.0069444444444444444444*h0*h1;      const double w54 = 0.0069444444444444444444*m_dx[0]*m_dx[1];
1877      const double w55 = 0.09672363354357992482*h0*h1;      const double w55 = 0.09672363354357992482*m_dx[0]*m_dx[1];
1878      const double w56 = 0.00049858867864229740201*h0*h1;      const double w56 = 0.00049858867864229740201*m_dx[0]*m_dx[1];
1879      const double w57 = 0.055555555555555555556*h0*h1;      const double w57 = 0.055555555555555555556*m_dx[0]*m_dx[1];
1880      const double w58 = 0.027777777777777777778*h0*h1;      const double w58 = 0.027777777777777777778*m_dx[0]*m_dx[1];
1881      const double w59 = 0.11111111111111111111*h0*h1;      const double w59 = 0.11111111111111111111*m_dx[0]*m_dx[1];
1882      const double w60 = -0.19716878364870322056*h1;      const double w60 = -0.19716878364870322056*m_dx[1];
1883      const double w61 = -0.19716878364870322056*h0;      const double w61 = -0.19716878364870322056*m_dx[0];
1884      const double w62 = -0.052831216351296779436*h0;      const double w62 = -0.052831216351296779436*m_dx[0];
1885      const double w63 = -0.052831216351296779436*h1;      const double w63 = -0.052831216351296779436*m_dx[1];
1886      const double w64 = 0.19716878364870322056*h1;      const double w64 = 0.19716878364870322056*m_dx[1];
1887      const double w65 = 0.052831216351296779436*h1;      const double w65 = 0.052831216351296779436*m_dx[1];
1888      const double w66 = 0.19716878364870322056*h0;      const double w66 = 0.19716878364870322056*m_dx[0];
1889      const double w67 = 0.052831216351296779436*h0;      const double w67 = 0.052831216351296779436*m_dx[0];
1890      const double w68 = -0.5*h1;      const double w68 = -0.5*m_dx[1];
1891      const double w69 = -0.5*h0;      const double w69 = -0.5*m_dx[0];
1892      const double w70 = 0.5*h1;      const double w70 = 0.5*m_dx[1];
1893      const double w71 = 0.5*h0;      const double w71 = 0.5*m_dx[0];
1894      const double w72 = 0.1555021169820365539*h0*h1;      const double w72 = 0.1555021169820365539*m_dx[0]*m_dx[1];
1895      const double w73 = 0.041666666666666666667*h0*h1;      const double w73 = 0.041666666666666666667*m_dx[0]*m_dx[1];
1896      const double w74 = 0.01116454968463011277*h0*h1;      const double w74 = 0.01116454968463011277*m_dx[0]*m_dx[1];
1897      const double w75 = 0.25*h0*h1;      const double w75 = 0.25*m_dx[0]*m_dx[1];
1898    
1899      rhs.requireWrite();      rhs.requireWrite();
1900  #pragma omp parallel  #pragma omp parallel
1901      {      {
1902          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1903  #pragma omp for  #pragma omp for
1904              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {              for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
1905                  for (index_t k0=0; k0<m_NE0; ++k0)  {                  for (index_t k0=0; k0<m_NE[0]; ++k0)  {
1906                      bool add_EM_S=false;                      bool add_EM_S=false;
1907                      bool add_EM_F=false;                      bool add_EM_F=false;
1908                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1909                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1910                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE[0]*k1;
1911                      ///////////////                      ///////////////
1912                      // process A //                      // process A //
1913                      ///////////////                      ///////////////
# Line 2162  void Rectangle::assemblePDESingle(Paso_S Line 2513  void Rectangle::assemblePDESingle(Paso_S
2513                      }                      }
2514    
2515                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2516                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
2517                      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);
2518                  } // end k0 loop                  } // end k0 loop
2519              } // end k1 loop              } // end k1 loop
# Line 2176  void Rectangle::assemblePDESingleReduced Line 2527  void Rectangle::assemblePDESingleReduced
2527          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2528          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
2529  {  {
2530      const double h0 = m_l0/m_gNE0;      const double w0 = -.25*m_dx[1]/m_dx[0];
     const double h1 = m_l1/m_gNE1;  
     const double w0 = -.25*h1/h0;  
2531      const double w1 = .25;      const double w1 = .25;
2532      const double w2 = -.25;      const double w2 = -.25;
2533      const double w3 = .25*h0/h1;      const double w3 = .25*m_dx[0]/m_dx[1];
2534      const double w4 = -.25*h0/h1;      const double w4 = -.25*m_dx[0]/m_dx[1];
2535      const double w5 = .25*h1/h0;      const double w5 = .25*m_dx[1]/m_dx[0];
2536      const double w6 = -.125*h1;      const double w6 = -.125*m_dx[1];
2537      const double w7 = -.125*h0;      const double w7 = -.125*m_dx[0];
2538      const double w8 = .125*h1;      const double w8 = .125*m_dx[1];
2539      const double w9 = .125*h0;      const double w9 = .125*m_dx[0];
2540      const double w10 = .0625*h0*h1;      const double w10 = .0625*m_dx[0]*m_dx[1];
2541      const double w11 = -.5*h1;      const double w11 = -.5*m_dx[1];
2542      const double w12 = -.5*h0;      const double w12 = -.5*m_dx[0];
2543      const double w13 = .5*h1;      const double w13 = .5*m_dx[1];
2544      const double w14 = .5*h0;      const double w14 = .5*m_dx[0];
2545      const double w15 = .25*h0*h1;      const double w15 = .25*m_dx[0]*m_dx[1];
2546    
2547      rhs.requireWrite();      rhs.requireWrite();
2548  #pragma omp parallel  #pragma omp parallel
2549      {      {
2550          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2551  #pragma omp for  #pragma omp for
2552              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {              for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
2553                  for (index_t k0=0; k0<m_NE0; ++k0)  {                  for (index_t k0=0; k0<m_NE[0]; ++k0)  {
2554                      bool add_EM_S=false;                      bool add_EM_S=false;
2555                      bool add_EM_F=false;                      bool add_EM_F=false;
2556                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
2557                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
2558                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE[0]*k1;
2559                      ///////////////                      ///////////////
2560                      // process A //                      // process A //
2561                      ///////////////                      ///////////////
# Line 2352  void Rectangle::assemblePDESingleReduced Line 2701  void Rectangle::assemblePDESingleReduced
2701                      }                      }
2702    
2703                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2704                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
2705                      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);
2706                  } // end k0 loop                  } // end k0 loop
2707              } // end k1 loop              } // end k1 loop
# Line 2366  void Rectangle::assemblePDESystem(Paso_S Line 2715  void Rectangle::assemblePDESystem(Paso_S
2715          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2716          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
2717  {  {
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
2718      dim_t numEq, numComp;      dim_t numEq, numComp;
2719      if (!mat)      if (!mat)
2720          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
# Line 2376  void Rectangle::assemblePDESystem(Paso_S Line 2723  void Rectangle::assemblePDESystem(Paso_S
2723          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
2724      }      }
2725    
2726      const double w0 = -0.1555021169820365539*h1/h0;      const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];
2727      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
2728      const double w2 = -0.15550211698203655390;      const double w2 = -0.15550211698203655390;
2729      const double w3 = 0.041666666666666666667*h0/h1;      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];
2730      const double w4 = 0.15550211698203655390;      const double w4 = 0.15550211698203655390;
2731      const double w5 = -0.041666666666666666667;      const double w5 = -0.041666666666666666667;
2732      const double w6 = -0.01116454968463011277*h1/h0;      const double w6 = -0.01116454968463011277*m_dx[1]/m_dx[0];
2733      const double w7 = 0.011164549684630112770;      const double w7 = 0.011164549684630112770;
2734      const double w8 = -0.011164549684630112770;      const double w8 = -0.011164549684630112770;
2735      const double w9 = -0.041666666666666666667*h1/h0;      const double w9 = -0.041666666666666666667*m_dx[1]/m_dx[0];
2736      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.041666666666666666667*m_dx[0]/m_dx[1];
2737      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = 0.1555021169820365539*m_dx[1]/m_dx[0];
2738      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = 0.1555021169820365539*m_dx[0]/m_dx[1];
2739      const double w13 = 0.01116454968463011277*h0/h1;      const double w13 = 0.01116454968463011277*m_dx[0]/m_dx[1];
2740      const double w14 = 0.01116454968463011277*h1/h0;      const double w14 = 0.01116454968463011277*m_dx[1]/m_dx[0];
2741      const double w15 = 0.041666666666666666667*h1/h0;      const double w15 = 0.041666666666666666667*m_dx[1]/m_dx[0];
2742      const double w16 = -0.01116454968463011277*h0/h1;      const double w16 = -0.01116454968463011277*m_dx[0]/m_dx[1];
2743      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*m_dx[0]/m_dx[1];
2744      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*m_dx[1]/m_dx[0];
2745      const double w19 = 0.25000000000000000000;      const double w19 = 0.25000000000000000000;
2746      const double w20 = -0.25000000000000000000;      const double w20 = -0.25000000000000000000;
2747      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*m_dx[0]/m_dx[1];
2748      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*m_dx[1]/m_dx[0];
2749      const double w23 = -0.16666666666666666667*h0/h1;      const double w23 = -0.16666666666666666667*m_dx[0]/m_dx[1];
2750      const double w24 = 0.33333333333333333333*h1/h0;      const double w24 = 0.33333333333333333333*m_dx[1]/m_dx[0];
2751      const double w25 = 0.33333333333333333333*h0/h1;      const double w25 = 0.33333333333333333333*m_dx[0]/m_dx[1];
2752      const double w26 = 0.16666666666666666667*h1/h0;      const double w26 = 0.16666666666666666667*m_dx[1]/m_dx[0];
2753      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = -0.33333333333333333333*m_dx[0]/m_dx[1];
2754      const double w28 = -0.032861463941450536761*h1;      const double w28 = -0.032861463941450536761*m_dx[1];
2755      const double w29 = -0.032861463941450536761*h0;      const double w29 = -0.032861463941450536761*m_dx[0];
2756      const double w30 = -0.12264065304058601714*h1;      const double w30 = -0.12264065304058601714*m_dx[1];
2757      const double w31 = -0.0023593469594139828636*h1;      const double w31 = -0.0023593469594139828636*m_dx[1];
2758      const double w32 = -0.008805202725216129906*h0;      const double w32 = -0.008805202725216129906*m_dx[0];
2759      const double w33 = -0.008805202725216129906*h1;      const double w33 = -0.008805202725216129906*m_dx[1];
2760      const double w34 = 0.032861463941450536761*h1;      const double w34 = 0.032861463941450536761*m_dx[1];
2761      const double w35 = 0.008805202725216129906*h1;      const double w35 = 0.008805202725216129906*m_dx[1];
2762      const double w36 = 0.008805202725216129906*h0;      const double w36 = 0.008805202725216129906*m_dx[0];
2763      const double w37 = 0.0023593469594139828636*h1;      const double w37 = 0.0023593469594139828636*m_dx[1];
2764      const double w38 = 0.12264065304058601714*h1;      const double w38 = 0.12264065304058601714*m_dx[1];
2765      const double w39 = 0.032861463941450536761*h0;      const double w39 = 0.032861463941450536761*m_dx[0];
2766      const double w40 = -0.12264065304058601714*h0;      const double w40 = -0.12264065304058601714*m_dx[0];
2767      const double w41 = -0.0023593469594139828636*h0;      const double w41 = -0.0023593469594139828636*m_dx[0];
2768      const double w42 = 0.0023593469594139828636*h0;      const double w42 = 0.0023593469594139828636*m_dx[0];
2769      const double w43 = 0.12264065304058601714*h0;      const double w43 = 0.12264065304058601714*m_dx[0];
2770      const double w44 = -0.16666666666666666667*h1;      const double w44 = -0.16666666666666666667*m_dx[1];
2771      const double w45 = -0.083333333333333333333*h0;      const double w45 = -0.083333333333333333333*m_dx[0];
2772      const double w46 = 0.083333333333333333333*h1;      const double w46 = 0.083333333333333333333*m_dx[1];
2773      const double w47 = 0.16666666666666666667*h1;      const double w47 = 0.16666666666666666667*m_dx[1];
2774      const double w48 = 0.083333333333333333333*h0;      const double w48 = 0.083333333333333333333*m_dx[0];
2775      const double w49 = -0.16666666666666666667*h0;      const double w49 = -0.16666666666666666667*m_dx[0];
2776      const double w50 = 0.16666666666666666667*h0;      const double w50 = 0.16666666666666666667*m_dx[0];
2777      const double w51 = -0.083333333333333333333*h1;      const double w51 = -0.083333333333333333333*m_dx[1];
2778      const double w52 = 0.025917019497006092316*h0*h1;      const double w52 = 0.025917019497006092316*m_dx[0]*m_dx[1];
2779      const double w53 = 0.0018607582807716854616*h0*h1;      const double w53 = 0.0018607582807716854616*m_dx[0]*m_dx[1];
2780      const double w54 = 0.0069444444444444444444*h0*h1;      const double w54 = 0.0069444444444444444444*m_dx[0]*m_dx[1];
2781      const double w55 = 0.09672363354357992482*h0*h1;      const double w55 = 0.09672363354357992482*m_dx[0]*m_dx[1];
2782      const double w56 = 0.00049858867864229740201*h0*h1;      const double w56 = 0.00049858867864229740201*m_dx[0]*m_dx[1];
2783      const double w57 = 0.055555555555555555556*h0*h1;      const double w57 = 0.055555555555555555556*m_dx[0]*m_dx[1];
2784      const double w58 = 0.027777777777777777778*h0*h1;      const double w58 = 0.027777777777777777778*m_dx[0]*m_dx[1];
2785      const double w59 = 0.11111111111111111111*h0*h1;      const double w59 = 0.11111111111111111111*m_dx[0]*m_dx[1];
2786      const double w60 = -0.19716878364870322056*h1;      const double w60 = -0.19716878364870322056*m_dx[1];
2787      const double w61 = -0.19716878364870322056*h0;      const double w61 = -0.19716878364870322056*m_dx[0];
2788      const double w62 = -0.052831216351296779436*h0;      const double w62 = -0.052831216351296779436*m_dx[0];
2789      const double w63 = -0.052831216351296779436*h1;      const double w63 = -0.052831216351296779436*m_dx[1];
2790      const double w64 = 0.19716878364870322056*h1;      const double w64 = 0.19716878364870322056*m_dx[1];
2791      const double w65 = 0.052831216351296779436*h1;      const double w65 = 0.052831216351296779436*m_dx[1];
2792      const double w66 = 0.19716878364870322056*h0;      const double w66 = 0.19716878364870322056*m_dx[0];
2793      const double w67 = 0.052831216351296779436*h0;      const double w67 = 0.052831216351296779436*m_dx[0];
2794      const double w68 = -0.5*h1;      const double w68 = -0.5*m_dx[1];
2795      const double w69 = -0.5*h0;      const double w69 = -0.5*m_dx[0];
2796      const double w70 = 0.5*h1;      const double w70 = 0.5*m_dx[1];
2797      const double w71 = 0.5*h0;      const double w71 = 0.5*m_dx[0];
2798      const double w72 = 0.1555021169820365539*h0*h1;      const double w72 = 0.1555021169820365539*m_dx[0]*m_dx[1];
2799      const double w73 = 0.041666666666666666667*h0*h1;      const double w73 = 0.041666666666666666667*m_dx[0]*m_dx[1];
2800      const double w74 = 0.01116454968463011277*h0*h1;      const double w74 = 0.01116454968463011277*m_dx[0]*m_dx[1];
2801      const double w75 = 0.25*h0*h1;      const double w75 = 0.25*m_dx[0]*m_dx[1];
2802    
2803      rhs.requireWrite();      rhs.requireWrite();
2804  #pragma omp parallel  #pragma omp parallel
2805      {      {
2806          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2807  #pragma omp for  #pragma omp for
2808              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {              for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
2809                  for (index_t k0=0; k0<m_NE0; ++k0)  {                  for (index_t k0=0; k0<m_NE[0]; ++k0)  {
2810                      bool add_EM_S=false;                      bool add_EM_S=false;
2811                      bool add_EM_F=false;                      bool add_EM_F=false;
2812                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
2813                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
2814                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE[0]*k1;
2815                      ///////////////                      ///////////////
2816                      // process A //                      // process A //
2817                      ///////////////                      ///////////////
# Line 3113  void Rectangle::assemblePDESystem(Paso_S Line 3460  void Rectangle::assemblePDESystem(Paso_S
3460                      }                      }
3461    
3462                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
3463                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
3464                      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,
3465                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
3466                  } // end k0 loop                  } // end k0 loop
# Line 3128  void Rectangle::assemblePDESystemReduced Line 3475  void Rectangle::assemblePDESystemReduced
3475          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
3476          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
3477  {  {
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
3478      dim_t numEq, numComp;      dim_t numEq, numComp;
3479      if (!mat)      if (!mat)
3480          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
# Line 3138  void Rectangle::assemblePDESystemReduced Line 3483  void Rectangle::assemblePDESystemReduced
3483          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
3484      }      }
3485    
3486      const double w0 = -.25*h1/h0;      const double w0 = -.25*m_dx[1]/m_dx[0];
3487      const double w1 = .25;      const double w1 = .25;
3488      const double w2 = -.25;      const double w2 = -.25;
3489      const double w3 = .25*h0/h1;      const double w3 = .25*m_dx[0]/m_dx[1];
3490      const double w4 = -.25*h0/h1;      const double w4 = -.25*m_dx[0]/m_dx[1];
3491      const double w5 = .25*h1/h0;      const double w5 = .25*m_dx[1]/m_dx[0];
3492      const double w6 = -.125*h1;      const double w6 = -.125*m_dx[1];
3493      const double w7 = -.125*h0;      const double w7 = -.125*m_dx[0];
3494      const double w8 = .125*h1;      const double w8 = .125*m_dx[1];
3495      const double w9 = .125*h0;      const double w9 = .125*m_dx[0];
3496      const double w10 = .0625*h0*h1;      const double w10 = .0625*m_dx[0]*m_dx[1];
3497      const double w11 = -.5*h1;      const double w11 = -.5*m_dx[1];
3498      const double w12 = -.5*h0;      const double w12 = -.5*m_dx[0];
3499      const double w13 = .5*h1;      const double w13 = .5*m_dx[1];
3500      const double w14 = .5*h0;      const double w14 = .5*m_dx[0];
3501      const double w15 = .25*h0*h1;      const double w15 = .25*m_dx[0]*m_dx[1];
3502    
3503      rhs.requireWrite();      rhs.requireWrite();
3504  #pragma omp parallel  #pragma omp parallel
3505      {      {
3506          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3507  #pragma omp for  #pragma omp for
3508              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {              for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
3509                  for (index_t k0=0; k0<m_NE0; ++k0)  {                  for (index_t k0=0; k0<m_NE[0]; ++k0)  {
3510                      bool add_EM_S=false;                      bool add_EM_S=false;
3511                      bool add_EM_F=false;                      bool add_EM_F=false;
3512                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
3513                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
3514                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE[0]*k1;
3515                      ///////////////                      ///////////////
3516                      // process A //                      // process A //
3517                      ///////////////                      ///////////////
# Line 3338  void Rectangle::assemblePDESystemReduced Line 3683  void Rectangle::assemblePDESystemReduced
3683                      }                      }
3684    
3685                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
3686                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
3687                      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,
3688                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
3689                  } // end k0 loop                  } // end k0 loop
# Line 3351  void Rectangle::assemblePDESystemReduced Line 3696  void Rectangle::assemblePDESystemReduced
3696  void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
3697        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3698  {  {
3699      const double h0 = m_l0/m_gNE0;      const double w0 = 0.31100423396407310779*m_dx[1];
3700      const double h1 = m_l1/m_gNE1;      const double w1 = 0.022329099369260225539*m_dx[1];
3701      const double w0 = 0.31100423396407310779*h1;      const double w10 = 0.022329099369260225539*m_dx[0];
3702      const double w1 = 0.022329099369260225539*h1;      const double w11 = 0.16666666666666666667*m_dx[0];
3703      const double w10 = 0.022329099369260225539*h0;      const double w12 = 0.33333333333333333333*m_dx[0];
3704      const double w11 = 0.16666666666666666667*h0;      const double w13 = 0.39433756729740644113*m_dx[0];
3705      const double w12 = 0.33333333333333333333*h0;      const double w14 = 0.10566243270259355887*m_dx[0];
3706      const double w13 = 0.39433756729740644113*h0;      const double w15 = 0.5*m_dx[0];
3707      const double w14 = 0.10566243270259355887*h0;      const double w2 = 0.083333333333333333333*m_dx[1];
3708      const double w15 = 0.5*h0;      const double w3 = 0.33333333333333333333*m_dx[1];
3709      const double w2 = 0.083333333333333333333*h1;      const double w4 = 0.16666666666666666667*m_dx[1];
3710      const double w3 = 0.33333333333333333333*h1;      const double w5 = 0.39433756729740644113*m_dx[1];
3711      const double w4 = 0.16666666666666666667*h1;      const double w6 = 0.10566243270259355887*m_dx[1];
3712      const double w5 = 0.39433756729740644113*h1;      const double w7 = 0.5*m_dx[1];
3713      const double w6 = 0.10566243270259355887*h1;      const double w8 = 0.083333333333333333333*m_dx[0];
3714      const double w7 = 0.5*h1;      const double w9 = 0.31100423396407310779*m_dx[0];
     const double w8 = 0.083333333333333333333*h0;  
     const double w9 = 0.31100423396407310779*h0;  
3715      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3716      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3717      rhs.requireWrite();      rhs.requireWrite();
# Line 3376  void Rectangle::assemblePDEBoundarySingl Line 3719  void Rectangle::assemblePDEBoundarySingl
3719      {      {
3720          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3721              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3722  #pragma omp for nowait  #pragma omp for
3723                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
3724                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3725                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3726                      const index_t e = k1;                      const index_t e = k1;
# Line 3430  void Rectangle::assemblePDEBoundarySingl Line 3773  void Rectangle::assemblePDEBoundarySingl
3773                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
3774                          }                          }
3775                      }                      }
3776                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_NN[0]*k1;
3777                      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);
3778                  }                  }
3779              } // end colouring              } // end colouring
# Line 3438  void Rectangle::assemblePDEBoundarySingl Line 3781  void Rectangle::assemblePDEBoundarySingl
3781    
3782          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
3783              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3784  #pragma omp for nowait  #pragma omp for
3785                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
3786                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3787                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3788                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
# Line 3492  void Rectangle::assemblePDEBoundarySingl Line 3835  void Rectangle::assemblePDEBoundarySingl
3835                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
3836                          }                          }
3837                      }                      }
3838                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_NN[0]*(k1+1)-2;
3839                      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);
3840                  }                  }
3841              } // end colouring              } // end colouring
# Line 3500  void Rectangle::assemblePDEBoundarySingl Line 3843  void Rectangle::assemblePDEBoundarySingl
3843    
3844          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
3845              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3846  #pragma omp for nowait  #pragma omp for
3847                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {
3848                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3849                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3850                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
# Line 3562  void Rectangle::assemblePDEBoundarySingl Line 3905  void Rectangle::assemblePDEBoundarySingl
3905    
3906          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
3907              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3908  #pragma omp for nowait  #pragma omp for
3909                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {
3910                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
3911                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3912                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3616  void Rectangle::assemblePDEBoundarySingl Line 3959  void Rectangle::assemblePDEBoundarySingl
3959                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
3960                          }                          }
3961                      }                      }
3962                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;
3963                      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);
3964                  }                  }
3965              } // end colouring              } // end colouring
# Line 3628  void Rectangle::assemblePDEBoundarySingl Line 3971  void Rectangle::assemblePDEBoundarySingl
3971  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
3972        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3973  {  {
3974      const double h0 = m_l0/m_gNE0;      const double w0 = 0.25*m_dx[1];
3975      const double h1 = m_l1/m_gNE1;      const double w1 = 0.5*m_dx[1];
3976      const double w0 = 0.25*h1;      const double w2 = 0.25*m_dx[0];
3977      const double w1 = 0.5*h1;      const double w3 = 0.5*m_dx[0];
     const double w2 = 0.25*h0;  
     const double w3 = 0.5*h0;  
3978      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3979      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3980      rhs.requireWrite();      rhs.requireWrite();
# Line 3641  void Rectangle::assemblePDEBoundarySingl Line 3982  void Rectangle::assemblePDEBoundarySingl
3982      {      {
3983          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3984              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3985  #pragma omp for nowait  #pragma omp for
3986                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
3987                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3988                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
3989                      const index_t e = k1;                      const index_t e = k1;
# Line 3667  void Rectangle::assemblePDEBoundarySingl Line 4008  void Rectangle::assemblePDEBoundarySingl
4008                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0_1;
4009                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
4010                      }                      }
4011                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_NN[0]*k1;
4012                      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);
4013                  }                  }
4014              } // end colouring              } // end colouring
# Line 3675  void Rectangle::assemblePDEBoundarySingl Line 4016  void Rectangle::assemblePDEBoundarySingl
4016    
4017          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
4018              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4019  #pragma omp for nowait  #pragma omp for
4020                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
4021                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
4022                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
4023                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
# Line 3701  void Rectangle::assemblePDEBoundarySingl Line 4042  void Rectangle::assemblePDEBoundarySingl
4042                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0_1;
4043                          EM_F[3]+=tmp0_1;                          EM_F[3]+=tmp0_1;
4044                      }                      }
4045                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_NN[0]*(k1+1)-2;
4046                      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);
4047                  }                  }
4048              } // end colouring              } // end colouring
# Line 3709  void Rectangle::assemblePDEBoundarySingl Line 4050  void Rectangle::assemblePDEBoundarySingl
4050    
4051          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
4052              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4053  #pragma omp for nowait  #pragma omp for
4054                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {
4055                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
4056                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
4057                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
# Line 3742  void Rectangle::assemblePDEBoundarySingl Line 4083  void Rectangle::assemblePDEBoundarySingl
4083    
4084          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4085              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4086  #pragma omp for nowait  #pragma omp for
4087                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {
4088                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
4089                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
4090                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
# Line 3767  void Rectangle::assemblePDEBoundarySingl Line 4108  void Rectangle::assemblePDEBoundarySingl
4108                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
4109                          EM_F[3]+=tmp0_1;                          EM_F[3]+=tmp0_1;
4110                      }                      }
4111                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;
4112                      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);
4113                  }                  }
4114              } // end colouring              } // end colouring
# Line 3779  void Rectangle::assemblePDEBoundarySingl Line 4120  void Rectangle::assemblePDEBoundarySingl
4120  void Rectangle::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
4121        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
4122  {  {
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
4123      dim_t numEq, numComp;      dim_t numEq, numComp;
4124      if (!mat) {      if (!mat) {
4125          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
# Line 3788  void Rectangle::assemblePDEBoundarySyste Line 4127  void Rectangle::assemblePDEBoundarySyste
4127          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
4128          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
4129      }      }
4130      const double w0 = 0.31100423396407310779*h1;      const double w0 = 0.31100423396407310779*m_dx[1];
4131      const double w1 = 0.022329099369260225539*h1;      const double w1 = 0.022329099369260225539*m_dx[1];
4132      const double w10 = 0.022329099369260225539*h0;      const double w10 = 0.022329099369260225539*m_dx[0];
4133      const double w11 = 0.16666666666666666667*h0;      const double w11 = 0.16666666666666666667*m_dx[0];
4134      const double w12 = 0.33333333333333333333*h0;      const double w12 = 0.33333333333333333333*m_dx[0];
4135      const double w13 = 0.39433756729740644113*h0;      const double w13 = 0.39433756729740644113*m_dx[0];
4136      const double w14 = 0.10566243270259355887*h0;      const double w14 = 0.10566243270259355887*m_dx[0];
4137      const double w15 = 0.5*h0;      const double w15 = 0.5*m_dx[0];
4138      const double w2 = 0.083333333333333333333*h1;      const double w2 = 0.083333333333333333333*m_dx[1];
4139      const double w3 = 0.33333333333333333333*h1;      const double w3 = 0.33333333333333333333*m_dx[1];
4140      const double w4 = 0.16666666666666666667*h1;      const double w4 = 0.16666666666666666667*m_dx[1];
4141      const double w5 = 0.39433756729740644113*h1;      const double w5 = 0.39433756729740644113*m_dx[1];
4142      const double w6 = 0.10566243270259355887*h1;      const double w6 = 0.10566243270259355887*m_dx[1];
4143      const double w7 = 0.5*h1;      const double w7 = 0.5*m_dx[1];
4144      const double w8 = 0.083333333333333333333*h0;      const double w8 = 0.083333333333333333333*m_dx[0];
4145      const double w9 = 0.31100423396407310779*h0;      const double w9 = 0.31100423396407310779*m_dx[0];
4146      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
4147      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
4148      rhs.requireWrite();      rhs.requireWrite();
# Line 3811  void Rectangle::assemblePDEBoundarySyste Line 4150  void Rectangle::assemblePDEBoundarySyste
4150      {      {
4151          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -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_NE[1]; 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 = k1;                      const index_t e = k1;
# Line 3877  void Rectangle::assemblePDEBoundarySyste Line 4216  void Rectangle::assemblePDEBoundarySyste
4216                              }                              }
4217                          }                          }
4218                      }                      }
4219                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_NN[0]*k1;
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);
4222                  }                  }
# Line 3886  void Rectangle::assemblePDEBoundarySyste Line 4225  void Rectangle::assemblePDEBoundarySyste
4225    
4226          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
4227              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4228  #pragma omp for nowait  #pragma omp for
4229                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=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[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
# Line 3952  void Rectangle::assemblePDEBoundarySyste Line 4291  void Rectangle::assemblePDEBoundarySyste
4291                              }                              }
4292                          }                          }
4293                      }                      }
4294                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_NN[0]*(k1+1)-2;
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);
4297                  }                  }
# Line 3961  void Rectangle::assemblePDEBoundarySyste Line 4300  void Rectangle::assemblePDEBoundarySyste
4300    
4301          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -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_NE[0]; 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[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
# Line 4036  void Rectangle::assemblePDEBoundarySyste Line 4375  void Rectangle::assemblePDEBoundarySyste
4375    
4376          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4377              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4378  #pragma omp for nowait  #pragma omp for
4379                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {
4380                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4381                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4382                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
# Line 4102  void Rectangle::assemblePDEBoundarySyste Line 4441  void Rectangle::assemblePDEBoundarySyste
4441                              }                              }
4442                          }                          }
4443                      }                      }
4444                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;
4445                      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,
4446                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
4447                  }                  }
# Line 4115  void Rectangle::assemblePDEBoundarySyste Line 4454  void Rectangle::assemblePDEBoundarySyste
4454  void Rectangle::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
4455        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
4456  {  {
     const double h0 = m_l0/m_gNE0;  
     const double h1 = m_l1/m_gNE1;  
4457      dim_t numEq, numComp;      dim_t numEq, numComp;
4458      if (!mat)      if (!mat)
4459          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());          numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
# Line 4124  void Rectangle::assemblePDEBoundarySyste Line 4461  void Rectangle::assemblePDEBoundarySyste
4461          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
4462          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
4463      }      }
4464      const double w0 = 0.25*h1;      const double w0 = 0.25*m_dx[1];
4465      const double w1 = 0.5*h1;      const double w1 = 0.5*m_dx[1];
4466      const double w2 = 0.25*h0;      const double w2 = 0.25*m_dx[0];
4467      const double w3 = 0.5*h0;      const double w3 = 0.5*m_dx[0];
4468      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
4469      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
4470      rhs.requireWrite();      rhs.requireWrite();
# Line 4135  void Rectangle::assemblePDEBoundarySyste Line 4472  void Rectangle::assemblePDEBoundarySyste
4472      {      {
4473          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
4474              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
4475  #pragma omp for nowait  #pragma omp for
4476                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
4477                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4478                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4479                      const index_t e = k1;                      const index_t e = k1;
# Line 4167  void Rectangle::assemblePDEBoundarySyste Line 4504  void Rectangle::assemblePDEBoundarySyste
4504                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
4505                          }                          }
4506                      }                      }
4507                      const index_t firstNode=m_N0*k1;                      const index_t firstNode=m_NN[0]*k1;
4508                      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,
4509                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
4510                  }                  }
# Line 4176  void Rectangle::assemblePDEBoundarySyste Line 4513  void Rectangle::assemblePDEBoundarySyste
4513    
4514          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
4515              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4516  #pragma omp for nowait  #pragma omp for
4517                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
4518                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4519                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4520                      const index_t e = m_faceOffset[1]+k1;                      const index_t e = m_faceOffset[1]+k1;
# Line 4208  void Rectangle::assemblePDEBoundarySyste Line 4545  void Rectangle::assemblePDEBoundarySyste
4545                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
4546                          }                          }
4547                      }                      }
4548                      const index_t firstNode=m_N0*(k1+1)-2;                      const index_t firstNode=m_NN[0]*(k1+1)-2;
4549                      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,
4550                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
4551                  }                  }
# Line 4217  void Rectangle::assemblePDEBoundarySyste Line 4554  void Rectangle::assemblePDEBoundarySyste
4554    
4555          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
4556              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4557  #pragma omp for nowait  #pragma omp for
4558                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {
4559                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4560                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4561                      const index_t e = m_faceOffset[2]+k0;                      const index_t e = m_faceOffset[2]+k0;
# Line 4258  void Rectangle::assemblePDEBoundarySyste Line 4595  void Rectangle::assemblePDEBoundarySyste
4595    
4596          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4597              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4598  #pragma omp for nowait  #pragma omp for
4599                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE[0]; k0+=2) {
4600                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4601                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
4602                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
# Line 4290  void Rectangle::assemblePDEBoundarySyste Line 4627  void Rectangle::assemblePDEBoundarySyste
4627                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
4628                          }                          }
4629                      }                      }
4630                      const index_t firstNode=m_N0*(m_N1-2)+k0;                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;
4631                      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,
4632                              firstNode, numEq, numComp);                              firstNode, numEq, numComp);
4633                  }                  }

Legend:
Removed from v.3793  
changed lines
  Added in v.4368

  ViewVC Help
Powered by ViewVC 1.1.26