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

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

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

branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp revision 3744 by caltinay, Tue Dec 13 06:41:54 2011 UTC trunk/ripley/src/Rectangle.cpp revision 4477 by caltinay, Thu Jun 20 02:41:28 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2011 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>
17  extern "C" {  #include <paso/SystemMatrix.h>
18  #include "paso/SystemMatrixPattern.h"  #include <esysUtils/esysFileWriter.h>
19  }  
20    #include <boost/scoped_array.hpp>
21    
22    #ifdef USE_NETCDF
23    #include <netcdfcpp.h>
24    #endif
25    
26  #if USE_SILO  #if USE_SILO
27  #include <silo.h>  #include <silo.h>
# Line 26  extern "C" { Line 33  extern "C" {
33  #include <iomanip>  #include <iomanip>
34    
35  using namespace std;  using namespace std;
36    using esysUtils::FileWriter;
37    
38  namespace ripley {  namespace ripley {
39    
40  Rectangle::Rectangle(int n0, int n1, double l0, double l1, int d0, int d1) :  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
41      RipleyDomain(2),                       double y1, int d0, int d1) :
42      m_gNE0(n0),      RipleyDomain(2)
43      m_gNE1(n1),  {
44      m_l0(l0),      // ignore subdivision parameters for serial run
45      m_l1(l1),      if (m_mpiInfo->size == 1) {
46      m_NX(d0),          d0=1;
47      m_NY(d1)          d1=1;
48  {      }
49    
50        bool warn=false;
51        // if number of subdivisions is non-positive, try to subdivide by the same
52        // ratio as the number of elements
53        if (d0<=0 && d1<=0) {
54            warn=true;
55            d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
56            d1=m_mpiInfo->size/d0;
57            if (d0*d1 != m_mpiInfo->size) {
58                // ratios not the same so subdivide side with more elements only
59                if (n0>n1) {
60                    d0=0;
61                    d1=1;
62                } else {
63                    d0=1;
64                    d1=0;
65                }
66            }
67        }
68        if (d0<=0) {
69            // d1 is preset, determine d0 - throw further down if result is no good
70            d0=m_mpiInfo->size/d1;
71        } else if (d1<=0) {
72            // d0 is preset, determine d1 - throw further down if result is no good
73            d1=m_mpiInfo->size/d0;
74        }
75    
76      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
77      // among number of ranks      // among number of ranks
78      if (m_NX*m_NY != m_mpiInfo->size)      if (d0*d1 != m_mpiInfo->size)
79          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
80    
81      if (n0%m_NX > 0 || n1%m_NY > 0)      if (warn) {
82          throw RipleyException("Number of elements must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
83                << d1 << "). This may not be optimal!" << endl;
84        }
85    
86        double l0 = x1-x0;
87        double l1 = y1-y0;
88        m_dx[0] = l0/n0;
89        m_dx[1] = l1/n1;
90    
91        if ((n0+1)%d0 > 0) {
92            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
93            l0=m_dx[0]*n0;
94            cout << "Warning: Adjusted number of elements and length. N0="
95                << n0 << ", l0=" << l0 << endl;
96        }
97        if ((n1+1)%d1 > 0) {
98            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
99            l1=m_dx[1]*n1;
100            cout << "Warning: Adjusted number of elements and length. N1="
101                << n1 << ", l1=" << l1 << endl;
102        }
103    
104        if ((d0 > 1 && (n0+1)/d0<2) || (d1 > 1 && (n1+1)/d1<2))
105            throw RipleyException("Too few elements for the number of ranks");
106    
107        m_gNE[0] = n0;
108        m_gNE[1] = n1;
109        m_origin[0] = x0;
110        m_origin[1] = y0;
111        m_length[0] = l0;
112        m_length[1] = l1;
113        m_NX[0] = d0;
114        m_NX[1] = d1;
115    
116        // local number of elements (with and without overlap)
117        m_NE[0] = m_ownNE[0] = (d0>1 ? (n0+1)/d0 : n0);
118        if (m_mpiInfo->rank%d0>0 && m_mpiInfo->rank%d0<d0-1)
119            m_NE[0]++;
120        else if (d0>1 && m_mpiInfo->rank%d0==d0-1)
121            m_ownNE[0]--;
122    
123        m_NE[1] = m_ownNE[1] = (d1>1 ? (n1+1)/d1 : n1);
124        if (m_mpiInfo->rank/d0>0 && m_mpiInfo->rank/d0<d1-1)
125            m_NE[1]++;
126        else if (d1>1 && m_mpiInfo->rank/d0==d1-1)
127            m_ownNE[1]--;
128    
129        // local number of nodes
130        m_NN[0] = m_NE[0]+1;
131        m_NN[1] = m_NE[1]+1;
132    
     // local number of elements  
     m_NE0 = n0/m_NX;  
     m_NE1 = n1/m_NY;  
     // local number of nodes (not necessarily owned)  
     m_N0 = m_NE0+1;  
     m_N1 = m_NE1+1;  
133      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
134      m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);      m_offset[0] = (n0+1)/d0*(m_mpiInfo->rank%d0);
135      m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);      if (m_offset[0] > 0)
136            m_offset[0]--;
137        m_offset[1] = (n1+1)/d1*(m_mpiInfo->rank/d0);
138        if (m_offset[1] > 0)
139            m_offset[1]--;
140    
141      populateSampleIds();      populateSampleIds();
142        createPattern();
143  }  }
144    
145  Rectangle::~Rectangle()  Rectangle::~Rectangle()
146  {  {
147        Paso_SystemMatrixPattern_free(m_pattern);
148        Paso_Connector_free(m_connector);
149  }  }
150    
151  string Rectangle::getDescription() const  string Rectangle::getDescription() const
# Line 72  bool Rectangle::operator==(const Abstrac Line 158  bool Rectangle::operator==(const Abstrac
158      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
159      if (o) {      if (o) {
160          return (RipleyDomain::operator==(other) &&          return (RipleyDomain::operator==(other) &&
161                  m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1                  m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1]
162                  && m_l0==o->m_l0 && m_l1==o->m_l1                  && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1]
163                  && m_NX==o->m_NX && m_NY==o->m_NY);                  && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1]
164                    && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1]);
165      }      }
166    
167      return false;      return false;
168  }  }
169    
170    void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
171                const vector<int>& first, const vector<int>& numValues,
172                const vector<int>& multiplier) const
173    {
174    #ifdef USE_NETCDF
175        // check destination function space
176        int myN0, myN1;
177        if (out.getFunctionSpace().getTypeCode() == Nodes) {
178            myN0 = m_NN[0];
179            myN1 = m_NN[1];
180        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
181                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
182            myN0 = m_NE[0];
183            myN1 = m_NE[1];
184        } else
185            throw RipleyException("readNcGrid(): invalid function space for output data object");
186    
187        if (first.size() != 2)
188            throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
189    
190        if (numValues.size() != 2)
191            throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
192    
193        if (multiplier.size() != 2)
194            throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");
195        for (size_t i=0; i<multiplier.size(); i++)
196            if (multiplier[i]<1)
197                throw RipleyException("readNcGrid(): all multipliers must be positive");
198    
199        // check file existence and size
200        NcFile f(filename.c_str(), NcFile::ReadOnly);
201        if (!f.is_valid())
202            throw RipleyException("readNcGrid(): cannot open file");
203    
204        NcVar* var = f.get_var(varname.c_str());
205        if (!var)
206            throw RipleyException("readNcGrid(): invalid variable");
207    
208        // TODO: rank>0 data support
209        const int numComp = out.getDataPointSize();
210        if (numComp > 1)
211            throw RipleyException("readNcGrid(): only scalar data supported");
212    
213        const int dims = var->num_dims();
214        boost::scoped_array<long> edges(var->edges());
215    
216        // is this a slice of the data object (dims!=2)?
217        // note the expected ordering of edges (as in numpy: y,x)
218        if ( (dims==2 && (numValues[1] > edges[0] || numValues[0] > edges[1]))
219                || (dims==1 && numValues[1]>1) ) {
220            throw RipleyException("readNcGrid(): not enough data in file");
221        }
222    
223        // check if this rank contributes anything
224        if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset[0] ||
225                first[1] >= m_offset[1]+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset[1])
226            return;
227    
228        // now determine how much this rank has to write
229    
230        // first coordinates in data object to write to
231        const int first0 = max(0, first[0]-m_offset[0]);
232        const int first1 = max(0, first[1]-m_offset[1]);
233        // indices to first value in file
234        const int idx0 = max(0, m_offset[0]-first[0]);
235        const int idx1 = max(0, m_offset[1]-first[1]);
236        // number of values to read
237        const int num0 = min(numValues[0]-idx0, myN0-first0);
238        const int num1 = min(numValues[1]-idx1, myN1-first1);
239    
240        vector<double> values(num0*num1);
241        if (dims==2) {
242            var->set_cur(idx1, idx0);
243            var->get(&values[0], num1, num0);
244        } else {
245            var->set_cur(idx0);
246            var->get(&values[0], num0);
247        }
248    
249        const int dpp = out.getNumDataPointsPerSample();
250        out.requireWrite();
251    
252        for (index_t y=0; y<num1; y++) {
253    #pragma omp parallel for
254            for (index_t x=0; x<num0; x++) {
255                const int baseIndex = first0+x*multiplier[0]
256                                      +(first1+y*multiplier[1])*myN0;
257                const int srcIndex = y*num0+x;
258                if (!isnan(values[srcIndex])) {
259                    for (index_t m1=0; m1<multiplier[1]; m1++) {
260                        for (index_t m0=0; m0<multiplier[0]; m0++) {
261                            const int dataIndex = baseIndex+m0+m1*myN0;
262                            double* dest = out.getSampleDataRW(dataIndex);
263                            for (index_t q=0; q<dpp; q++) {
264                                *dest++ = values[srcIndex];
265                            }
266                        }
267                    }
268                }
269            }
270        }
271    #else
272        throw RipleyException("readNcGrid(): not compiled with netCDF support");
273    #endif
274    }
275    
276    void Rectangle::readBinaryGrid(escript::Data& out, string filename,
277                                   const vector<int>& first,
278                                   const vector<int>& numValues,
279                                   const vector<int>& multiplier) const
280    {
281        // check destination function space
282        int myN0, myN1;
283        if (out.getFunctionSpace().getTypeCode() == Nodes) {
284            myN0 = m_NN[0];
285            myN1 = m_NN[1];
286        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
287                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
288            myN0 = m_NE[0];
289            myN1 = m_NE[1];
290        } else
291            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
292    
293        // check file existence and size
294        ifstream f(filename.c_str(), ifstream::binary);
295        if (f.fail()) {
296            throw RipleyException("readBinaryGrid(): cannot open file");
297        }
298        f.seekg(0, ios::end);
299        const int numComp = out.getDataPointSize();
300        const int filesize = f.tellg();
301        const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(float);
302        if (filesize < reqsize) {
303            f.close();
304            throw RipleyException("readBinaryGrid(): not enough data in file");
305        }
306    
307        // check if this rank contributes anything
308        if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0] <= m_offset[0] ||
309                first[1] >= m_offset[1]+myN1 || first[1]+numValues[1] <= m_offset[1]) {
310            f.close();
311            return;
312        }
313    
314        // now determine how much this rank has to write
315    
316        // first coordinates in data object to write to
317        const int first0 = max(0, first[0]-m_offset[0]);
318        const int first1 = max(0, first[1]-m_offset[1]);
319        // indices to first value in file
320        const int idx0 = max(0, m_offset[0]-first[0]);
321        const int idx1 = max(0, m_offset[1]-first[1]);
322        // number of values to read
323        const int num0 = min(numValues[0]-idx0, myN0-first0);
324        const int num1 = min(numValues[1]-idx1, myN1-first1);
325    
326        out.requireWrite();
327        vector<float> values(num0*numComp);
328        const int dpp = out.getNumDataPointsPerSample();
329    
330        for (index_t y=0; y<num1; y++) {
331            const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);
332            f.seekg(fileofs*sizeof(float));
333            f.read((char*)&values[0], num0*numComp*sizeof(float));
334            for (index_t x=0; x<num0; x++) {
335                const int baseIndex = first0+x*multiplier[0]
336                                        +(first1+y*multiplier[1])*myN0;
337                for (index_t m1=0; m1<multiplier[1]; m1++) {
338                    for (index_t m0=0; m0<multiplier[0]; m0++) {
339                        const int dataIndex = baseIndex+m0+m1*myN0;
340                        double* dest = out.getSampleDataRW(dataIndex);
341                        for (index_t c=0; c<numComp; c++) {
342                            if (!std::isnan(values[x*numComp+c])) {
343                                for (index_t q=0; q<dpp; q++) {
344                                    *dest++ = static_cast<double>(values[x*numComp+c]);
345                                }
346                            }
347                        }
348                    }
349                }
350            }
351        }
352    
353        f.close();
354    }
355    
356    void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,
357                                    int byteOrder, int dataType) const
358    {
359        // the mapping is not universally correct but should work on our
360        // supported platforms
361        switch (dataType) {
362            case DATATYPE_INT32:
363                writeBinaryGridImpl<int>(in, filename, byteOrder);
364                break;
365            case DATATYPE_FLOAT32:
366                writeBinaryGridImpl<float>(in, filename, byteOrder);
367                break;
368            case DATATYPE_FLOAT64:
369                writeBinaryGridImpl<double>(in, filename, byteOrder);
370                break;
371            default:
372                throw RipleyException("writeBinaryGrid(): invalid or unsupported datatype");
373        }
374    }
375    
376    template<typename ValueType>
377    void Rectangle::writeBinaryGridImpl(const escript::Data& in,
378                                        const string& filename, int byteOrder) const
379    {
380        // check function space and determine number of points
381        int myN0, myN1;
382        int totalN0, totalN1;
383        if (in.getFunctionSpace().getTypeCode() == Nodes) {
384            myN0 = m_NN[0];
385            myN1 = m_NN[1];
386            totalN0 = m_gNE[0]+1;
387            totalN1 = m_gNE[1]+1;
388        } else if (in.getFunctionSpace().getTypeCode() == Elements ||
389                    in.getFunctionSpace().getTypeCode() == ReducedElements) {
390            myN0 = m_NE[0];
391            myN1 = m_NE[1];
392            totalN0 = m_gNE[0];
393            totalN1 = m_gNE[1];
394        } else
395            throw RipleyException("writeBinaryGrid(): invalid function space of data object");
396    
397        const int numComp = in.getDataPointSize();
398        const int dpp = in.getNumDataPointsPerSample();
399    
400        if (numComp > 1 || dpp > 1)
401            throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
402    
403        escript::Data* _in = const_cast<escript::Data*>(&in);
404        const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;
405    
406        // from here on we know that each sample consists of one value
407        FileWriter* fw = new FileWriter();
408        fw->openFile(filename, fileSize);
409        MPIBarrier();
410    
411        for (index_t y=0; y<myN1; y++) {
412            const int fileofs = (m_offset[0]+(m_offset[1]+y)*totalN0)*sizeof(ValueType);
413            ostringstream oss;
414    
415            for (index_t x=0; x<myN0; x++) {
416                const double* sample = _in->getSampleDataRO(y*myN0+x);
417                ValueType fvalue = static_cast<ValueType>(*sample);
418                if (byteOrder == BYTEORDER_NATIVE) {
419                    oss.write((char*)&fvalue, sizeof(fvalue));
420                } else {
421                    char* value = reinterpret_cast<char*>(&fvalue);
422                    oss.write(byte_swap32(value), sizeof(fvalue));
423                }
424            }
425            fw->writeAt(oss, fileofs);
426        }
427        fw->close();
428    }
429    
430  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
431  {  {
432  #if USE_SILO  #if USE_SILO
# Line 88  void Rectangle::dump(const string& fileN Line 435  void Rectangle::dump(const string& fileN
435          fn+=".silo";          fn+=".silo";
436      }      }
437    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
438      int driver=DB_HDF5;          int driver=DB_HDF5;    
     string siloPath;  
439      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
440        const char* blockDirFmt = "/block%04d";
441    
442  #ifdef ESYS_MPI  #ifdef ESYS_MPI
443      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
444        const int NUM_SILO_FILES = 1;
445  #endif  #endif
446    
447      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 111  void Rectangle::dump(const string& fileN Line 457  void Rectangle::dump(const string& fileN
457                          PMPIO_DefaultClose, (void*)&driver);                          PMPIO_DefaultClose, (void*)&driver);
458          }          }
459          if (baton) {          if (baton) {
460              char str[64];              char siloPath[64];
461              snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));              snprintf(siloPath, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
462              siloPath = str;              dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath);
             dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());  
463          }          }
464  #endif  #endif
465      } else {      } else {
# Line 126  void Rectangle::dump(const string& fileN Line 471  void Rectangle::dump(const string& fileN
471              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
472                      getDescription().c_str(), driver);                      getDescription().c_str(), driver);
473          }          }
474            char siloPath[64];
475            snprintf(siloPath, 64, blockDirFmt, 0);
476            DBMkDir(dbfile, siloPath);
477            DBSetDir(dbfile, siloPath);
478      }      }
479    
480      if (!dbfile)      if (!dbfile)
# Line 140  void Rectangle::dump(const string& fileN Line 489  void Rectangle::dump(const string& fileN
489      }      }
490      */      */
491    
492      boost::scoped_ptr<double> x(new double[m_N0]);      boost::scoped_ptr<double> x(new double[m_NN[0]]);
493      boost::scoped_ptr<double> y(new double[m_N1]);      boost::scoped_ptr<double> y(new double[m_NN[1]]);
494      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);  
495  #pragma omp parallel  #pragma omp parallel
496      {      {
497  #pragma omp for nowait  #pragma omp for nowait
498          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
499              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=getLocalCoordinate(i0, 0);
500          }          }
501  #pragma omp for nowait  #pragma omp for nowait
502          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
503              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=getLocalCoordinate(i1, 1);
504          }          }
505      }      }
506      IndexVector dims = getNumNodesPerDim();      int* dims = const_cast<int*>(getNumNodesPerDim());
507    
508      // write mesh      // write mesh
509      DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,
510              DB_COLLINEAR, NULL);              DB_COLLINEAR, NULL);
511    
512      // write node ids      // write node ids
513      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 2,      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 2,
514              NULL, 0, DB_INT, DB_NODECENT, NULL);              NULL, 0, DB_INT, DB_NODECENT, NULL);
515    
516      // write element ids      // write element ids
517      dims = getNumElementsPerDim();      dims = const_cast<int*>(getNumElementsPerDim());
518      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
519              &dims[0], 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);              dims, 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);
520    
521      // rank 0 writes multimesh and multivar      // rank 0 writes multimesh and multivar
522      if (m_mpiInfo->rank == 0) {      if (m_mpiInfo->rank == 0) {
# Line 218  void Rectangle::dump(const string& fileN Line 565  void Rectangle::dump(const string& fileN
565      }      }
566    
567  #else // USE_SILO  #else // USE_SILO
568      throw RipleyException("dump(): no Silo support");      throw RipleyException("dump: no Silo support");
569  #endif  #endif
570  }  }
571    
# Line 226  const int* Rectangle::borrowSampleRefere Line 573  const int* Rectangle::borrowSampleRefere
573  {  {
574      switch (fsType) {      switch (fsType) {
575          case Nodes:          case Nodes:
576            case ReducedNodes: // FIXME: reduced
577              return &m_nodeId[0];              return &m_nodeId[0];
578          case DegreesOfFreedom:          case DegreesOfFreedom:
579          case ReducedDegreesOfFreedom: //FIXME: reduced          case ReducedDegreesOfFreedom: // FIXME: reduced
580              return &m_dofId[0];              return &m_dofId[0];
581          case Elements:          case Elements:
582          case ReducedElements:          case ReducedElements:
# Line 241  const int* Rectangle::borrowSampleRefere Line 589  const int* Rectangle::borrowSampleRefere
589      }      }
590    
591      stringstream msg;      stringstream msg;
592      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs: invalid function space type " << fsType;
         << functionSpaceTypeAsString(fsType);  
593      throw RipleyException(msg.str());      throw RipleyException(msg.str());
594  }  }
595    
596  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsType, index_t id) const
597  {  {
598  #ifdef ESYS_MPI      if (getMPISize()==1)
599      if (fsCode != ReducedNodes) {          return true;
600          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];  
601          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;      switch (fsType) {
602          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);          case Nodes:
603            case ReducedNodes: // FIXME: reduced
604                return (m_dofMap[id] < getNumDOF());
605            case DegreesOfFreedom:
606            case ReducedDegreesOfFreedom:
607                return true;
608            case Elements:
609            case ReducedElements:
610                // check ownership of element's bottom left node
611                return (m_dofMap[id%m_NE[0]+m_NN[0]*(id/m_NE[0])] < getNumDOF());
612            case FaceElements:
613            case ReducedFaceElements:
614                {
615                    // determine which face the sample belongs to before
616                    // checking ownership of corresponding element's first node
617                    dim_t n=0;
618                    for (size_t i=0; i<4; i++) {
619                        n+=m_faceCount[i];
620                        if (id<n) {
621                            index_t k;
622                            if (i==1)
623                                k=m_NN[0]-2;
624                            else if (i==3)
625                                k=m_NN[0]*(m_NN[1]-2);
626                            else
627                                k=0;
628                            // determine whether to move right or up
629                            const index_t delta=(i/2==0 ? m_NN[0] : 1);
630                            return (m_dofMap[k+(id-n+m_faceCount[i])*delta] < getNumDOF());
631                        }
632                    }
633                    return false;
634                }
635            default:
636                break;
637        }
638    
639        stringstream msg;
640        msg << "ownSample: invalid function space type " << fsType;
641        throw RipleyException(msg.str());
642    }
643    
644    void Rectangle::setToNormal(escript::Data& out) const
645    {
646        if (out.getFunctionSpace().getTypeCode() == FaceElements) {
647            out.requireWrite();
648    #pragma omp parallel
649            {
650                if (m_faceOffset[0] > -1) {
651    #pragma omp for nowait
652                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
653                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
654                        // set vector at two quadrature points
655                        *o++ = -1.;
656                        *o++ = 0.;
657                        *o++ = -1.;
658                        *o = 0.;
659                    }
660                }
661    
662                if (m_faceOffset[1] > -1) {
663    #pragma omp for nowait
664                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
665                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
666                        // set vector at two quadrature points
667                        *o++ = 1.;
668                        *o++ = 0.;
669                        *o++ = 1.;
670                        *o = 0.;
671                    }
672                }
673    
674                if (m_faceOffset[2] > -1) {
675    #pragma omp for nowait
676                    for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
677                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
678                        // set vector at two quadrature points
679                        *o++ = 0.;
680                        *o++ = -1.;
681                        *o++ = 0.;
682                        *o = -1.;
683                    }
684                }
685    
686                if (m_faceOffset[3] > -1) {
687    #pragma omp for nowait
688                    for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
689                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
690                        // set vector at two quadrature points
691                        *o++ = 0.;
692                        *o++ = 1.;
693                        *o++ = 0.;
694                        *o = 1.;
695                    }
696                }
697            } // end of parallel section
698        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
699            out.requireWrite();
700    #pragma omp parallel
701            {
702                if (m_faceOffset[0] > -1) {
703    #pragma omp for nowait
704                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
705                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
706                        *o++ = -1.;
707                        *o = 0.;
708                    }
709                }
710    
711                if (m_faceOffset[1] > -1) {
712    #pragma omp for nowait
713                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
714                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
715                        *o++ = 1.;
716                        *o = 0.;
717                    }
718                }
719    
720                if (m_faceOffset[2] > -1) {
721    #pragma omp for nowait
722                    for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
723                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
724                        *o++ = 0.;
725                        *o = -1.;
726                    }
727                }
728    
729                if (m_faceOffset[3] > -1) {
730    #pragma omp for nowait
731                    for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
732                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
733                        *o++ = 0.;
734                        *o = 1.;
735                    }
736                }
737            } // end of parallel section
738    
739      } else {      } else {
740          stringstream msg;          stringstream msg;
741          msg << "ownSample() not implemented for "          msg << "setToNormal: invalid function space type "
742              << functionSpaceTypeAsString(fsCode);              << out.getFunctionSpace().getTypeCode();
743          throw RipleyException(msg.str());          throw RipleyException(msg.str());
744      }      }
 #else  
     return true;  
 #endif  
745  }  }
746    
747  void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const  void Rectangle::setToSize(escript::Data& out) const
748    {
749        if (out.getFunctionSpace().getTypeCode() == Elements
750                || out.getFunctionSpace().getTypeCode() == ReducedElements) {
751            out.requireWrite();
752            const dim_t numQuad=out.getNumDataPointsPerSample();
753            const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]);
754    #pragma omp parallel for
755            for (index_t k = 0; k < getNumElements(); ++k) {
756                double* o = out.getSampleDataRW(k);
757                fill(o, o+numQuad, size);
758            }
759        } else if (out.getFunctionSpace().getTypeCode() == FaceElements
760                || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
761            out.requireWrite();
762            const dim_t numQuad=out.getNumDataPointsPerSample();
763    #pragma omp parallel
764            {
765                if (m_faceOffset[0] > -1) {
766    #pragma omp for nowait
767                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
768                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
769                        fill(o, o+numQuad, m_dx[1]);
770                    }
771                }
772    
773                if (m_faceOffset[1] > -1) {
774    #pragma omp for nowait
775                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
776                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
777                        fill(o, o+numQuad, m_dx[1]);
778                    }
779                }
780    
781                if (m_faceOffset[2] > -1) {
782    #pragma omp for nowait
783                    for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
784                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
785                        fill(o, o+numQuad, m_dx[0]);
786                    }
787                }
788    
789                if (m_faceOffset[3] > -1) {
790    #pragma omp for nowait
791                    for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
792                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
793                        fill(o, o+numQuad, m_dx[0]);
794                    }
795                }
796            } // end of parallel section
797    
798        } else {
799            stringstream msg;
800            msg << "setToSize: invalid function space type "
801                << out.getFunctionSpace().getTypeCode();
802            throw RipleyException(msg.str());
803        }
804    }
805    
806    void Rectangle::Print_Mesh_Info(const bool full) const
807    {
808        RipleyDomain::Print_Mesh_Info(full);
809        if (full) {
810            cout << "     Id  Coordinates" << endl;
811            cout.precision(15);
812            cout.setf(ios::scientific, ios::floatfield);
813            for (index_t i=0; i < getNumNodes(); i++) {
814                cout << "  " << setw(5) << m_nodeId[i]
815                    << "  " << getLocalCoordinate(i%m_NN[0], 0)
816                    << "  " << getLocalCoordinate(i/m_NN[0], 1) << endl;
817            }
818        }
819    }
820    
821    
822    //protected
823    void Rectangle::assembleCoordinates(escript::Data& arg) const
824    {
825        escriptDataC x = arg.getDataC();
826        int numDim = m_numDim;
827        if (!isDataPointShapeEqual(&x, 1, &numDim))
828            throw RipleyException("setToX: Invalid Data object shape");
829        if (!numSamplesEqual(&x, 1, getNumNodes()))
830            throw RipleyException("setToX: Illegal number of samples in Data object");
831    
832        arg.requireWrite();
833    #pragma omp parallel for
834        for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
835            for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
836                double* point = arg.getSampleDataRW(i0+m_NN[0]*i1);
837                point[0] = getLocalCoordinate(i0, 0);
838                point[1] = getLocalCoordinate(i1, 1);
839            }
840        }
841    }
842    
843    //protected
844    void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
845  {  {
     escript::Data& in = *const_cast<escript::Data*>(&cIn);  
846      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
847      const double h0 = m_l0/m_gNE0;      const double cx0 = .21132486540518711775/m_dx[0];
848      const double h1 = m_l1/m_gNE1;      const double cx1 = .78867513459481288225/m_dx[0];
849      const double cx0 = -1./h0;      const double cx2 = 1./m_dx[0];
850      const double cx1 = -.78867513459481288225/h0;      const double cy0 = .21132486540518711775/m_dx[1];
851      const double cx2 = -.5/h0;      const double cy1 = .78867513459481288225/m_dx[1];
852      const double cx3 = -.21132486540518711775/h0;      const double cy2 = 1./m_dx[1];
     const double cx4 = .21132486540518711775/h0;  
     const double cx5 = .5/h0;  
     const double cx6 = .78867513459481288225/h0;  
     const double cx7 = 1./h0;  
     const double cy0 = -1./h1;  
     const double cy1 = -.78867513459481288225/h1;  
     const double cy2 = -.5/h1;  
     const double cy3 = -.21132486540518711775/h1;  
     const double cy4 = .21132486540518711775/h1;  
     const double cy5 = .5/h1;  
     const double cy6 = .78867513459481288225/h1;  
     const double cy7 = 1./h1;  
853    
854      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
855          /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */          out.requireWrite();
856  #pragma omp parallel for  #pragma omp parallel
857          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
858              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
859                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
860                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
861                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
862                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
863                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
864                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
865                      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));
866                      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));
867                      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));
868                      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));
869                      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]));
870                      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) {
871                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
872                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
873                  } /* end of component loop i */                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
874              } /* end of k0 loop */                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
875          } /* end of k1 loop */                          o[INDEX3(i,0,2,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
876          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */                          o[INDEX3(i,1,2,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
877                            o[INDEX3(i,0,3,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
878                            o[INDEX3(i,1,3,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
879                        } // end of component loop i
880                    } // end of k0 loop
881                } // end of k1 loop
882            } // end of parallel section
883      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
884          /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */          out.requireWrite();
885  #pragma omp parallel for  #pragma omp parallel
886          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
887              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
888                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
889                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
890                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
891                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
892                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
893                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
894                      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));
895                      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));
896                  } /* end of component loop i */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
897              } /* end of k0 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
898          } /* end of k1 loop */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
899          /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */                      for (index_t i=0; i < numComp; ++i) {
900                            o[INDEX3(i,0,0,numComp,2)] = (f_10[i] + f_11[i] - f_00[i] - f_01[i])*cx2/2;
901                            o[INDEX3(i,1,0,numComp,2)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
902                        } // end of component loop i
903                    } // end of k0 loop
904                } // end of k1 loop
905            } // end of parallel section
906      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
907            out.requireWrite();
908  #pragma omp parallel  #pragma omp parallel
909          {          {
910              /*** GENERATOR SNIP_GRAD_FACES TOP */              vector<double> f_00(numComp);
911                vector<double> f_01(numComp);
912                vector<double> f_10(numComp);
913                vector<double> f_11(numComp);
914              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
915  #pragma omp for nowait  #pragma omp for nowait
916                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
917                      const register 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));
918                      const register 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));
919                      const register 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));
920                      const register 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));
921                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
922                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
923                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
924                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy2;
925                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
926                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy2;
927                      } /* end of component loop i */                      } // end of component loop i
928                  } /* end of k1 loop */                  } // end of k1 loop
929              } /* end of face 0 */              } // end of face 0
930              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
931  #pragma omp for nowait  #pragma omp for nowait
932                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
933                      const register 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));
934                      const register 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));
935                      const register 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));
936                      const register 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));
937                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
938                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
939                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx1 + (f_11[i]-f_01[i])*cx0;
940                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_11[i]-f_10[i])*cy2;
941                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx0 + (f_11[i]-f_01[i])*cx1;
942                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = (f_11[i]-f_10[i])*cy2;
943                      } /* end of component loop i */                      } // end of component loop i
944                  } /* end of k1 loop */                  } // end of k1 loop
945              } /* end of face 1 */              } // end of face 1
946              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
947  #pragma omp for nowait  #pragma omp for nowait
948                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
949                      const register 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));
950                      const register 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));
951                      const register 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));
952                      const register 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));
953                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
954                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
955                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx2;
956                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
957                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = (f_10[i]-f_00[i])*cx2;
958                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
959                      } /* end of component loop i */                      } // end of component loop i
960                  } /* end of k0 loop */                  } // end of k0 loop
961              } /* end of face 2 */              } // end of face 2
962              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
963  #pragma omp for nowait  #pragma omp for nowait
964                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
965                      const register 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));
966                      const register 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));
967                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
968                      const register 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));
969                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
970                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
971                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_11[i]-f_01[i])*cx2;
972                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy1 + (f_11[i]-f_10[i])*cy0;
973                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = (f_11[i]-f_01[i])*cx2;
974                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = (f_01[i]-f_00[i])*cy0 + (f_11[i]-f_10[i])*cy1;
975                      } /* end of component loop i */                      } // end of component loop i
976                  } /* end of k0 loop */                  } // end of k0 loop
977              } /* end of face 3 */              } // end of face 3
             /* GENERATOR SNIP_GRAD_FACES BOTTOM */  
978          } // end of parallel section          } // end of parallel section
979    
980      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
981            out.requireWrite();
982  #pragma omp parallel  #pragma omp parallel
983          {          {
984              /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */              vector<double> f_00(numComp);
985                vector<double> f_01(numComp);
986                vector<double> f_10(numComp);
987                vector<double> f_11(numComp);
988              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
989  #pragma omp for nowait  #pragma omp for nowait
990                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
991                      const register 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));
992                      const register 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));
993                      const register 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));
994                      const register 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));
995                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
996                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
997                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i] + f_11[i] - f_00[i] - f_01[i])*cx2/2;
998                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i]-f_00[i])*cy2;
999                      } /* end of component loop i */                      } // end of component loop i
1000                  } /* end of k1 loop */                  } // end of k1 loop
1001              } /* end of face 0 */              } // end of face 0
1002              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1003  #pragma omp for nowait  #pragma omp for nowait
1004                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1005                      const register 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));
1006                      const register 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));
1007                      const register 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));
1008                      const register 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));
1009                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1010                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1011                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i] + f_11[i] - f_00[i] - f_01[i])*cx2/2;
1012                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = (f_11[i]-f_10[i])*cy2;
1013                      } /* end of component loop i */                      } // end of component loop i
1014                  } /* end of k1 loop */                  } // end of k1 loop
1015              } /* end of face 1 */              } // end of face 1
1016              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1017  #pragma omp for nowait  #pragma omp for nowait
1018                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1019                      const register 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));
1020                      const register 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));
1021                      const register 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));
1022                      const register 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));
1023                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1024                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1025                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_10[i]-f_00[i])*cx2;
1026                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
1027                      } /* end of component loop i */                      } // end of component loop i
1028                  } /* end of k0 loop */                  } // end of k0 loop
1029              } /* end of face 2 */              } // end of face 2
1030              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1031  #pragma omp for nowait  #pragma omp for nowait
1032                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1033                      const register 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));
1034                      const register 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));
1035                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-2, m_NN[0])), numComp*sizeof(double));
1036                      const register 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));
1037                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1038                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1039                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = (f_11[i]-f_01[i])*cx2;
1040                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);                          o[INDEX3(i,1,0,numComp,2)] = (f_01[i] + f_11[i] - f_00[i] - f_10[i])*cy2/2;
1041                      } /* end of component loop i */                      } // end of component loop i
1042                  } /* end of k0 loop */                  } // end of k0 loop
1043              } /* end of face 3 */              } // end of face 3
             /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */  
1044          } // end of parallel section          } // end of parallel section
     } else {  
         stringstream msg;  
         msg << "setToGradient() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
1045      }      }
1046  }  }
1047    
1048  void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  //protected
1049    void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1050  {  {
1051      escript::Data& in = *const_cast<escript::Data*>(&arg);      const dim_t numComp = arg.getDataPointSize();
1052      const dim_t numComp = in.getDataPointSize();      const index_t left = (m_offset[0]==0 ? 0 : 1);
1053      const double h0 = m_l0/m_gNE0;      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1054      const double h1 = m_l1/m_gNE1;      const int fs=arg.getFunctionSpace().getTypeCode();
1055      if (arg.getFunctionSpace().getTypeCode() == Elements) {      if (fs == Elements && arg.actsExpanded()) {
         const double w_0 = h0*h1/4.;  
1056  #pragma omp parallel  #pragma omp parallel
1057          {          {
1058              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1059                const double w = m_dx[0]*m_dx[1]/4.;
1060  #pragma omp for nowait  #pragma omp for nowait
1061              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1062                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1063                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1064                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1065                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1066                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1067                          const register double f_2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
1068                          const register double f_3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
1069                          int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;                          int_local[i]+=(f0+f1+f2+f3)*w;
1070                      }  /* end of component loop i */                      }  // end of component loop i
1071                  } /* end of k0 loop */                  } // end of k0 loop
1072              } /* end of k1 loop */              } // end of k1 loop
   
1073  #pragma omp critical  #pragma omp critical
1074              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1075                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1076          } // end of parallel section          } // end of parallel section
1077      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
1078          const double w_0 = h0*h1;      } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1079            const double w = m_dx[0]*m_dx[1];
1080  #pragma omp parallel  #pragma omp parallel
1081          {          {
1082              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1083  #pragma omp for nowait  #pragma omp for nowait
1084              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1085                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1086                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE[0]));
1087                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1088                          int_local[i]+=f[i]*w_0;                          int_local[i]+=f[i]*w;
1089                      }  /* end of component loop i */                      }
1090                  } /* end of k0 loop */                  }
1091              } /* end of k1 loop */              }
   
1092  #pragma omp critical  #pragma omp critical
1093              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1094                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1095          } // end of parallel section          } // end of parallel section
1096      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
1097          const double w_0 = h0/2.;      } else if (fs == FaceElements && arg.actsExpanded()) {
         const double w_1 = h1/2.;  
1098  #pragma omp parallel  #pragma omp parallel
1099          {          {
1100              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1101                const double w0 = m_dx[0]/2.;
1102                const double w1 = m_dx[1]/2.;
1103              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1104  #pragma omp for nowait  #pragma omp for nowait
1105                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1106                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1107                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1108                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1109                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1110                          int_local[i]+=(f_0+f_1)*w_1;                          int_local[i]+=(f0+f1)*w1;
1111                      }  /* end of component loop i */                      }  // end of component loop i
1112                  } /* end of k1 loop */                  } // end of k1 loop
1113              }              }
1114    
1115              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1116  #pragma omp for nowait  #pragma omp for nowait
1117                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1118                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1119                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1120                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1121                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1122                          int_local[i]+=(f_0+f_1)*w_1;                          int_local[i]+=(f0+f1)*w1;
1123                      }  /* end of component loop i */                      }  // end of component loop i
1124                  } /* end of k1 loop */                  } // end of k1 loop
1125              }              }
1126    
1127              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1128  #pragma omp for nowait  #pragma omp for nowait
1129                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1130                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1131                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1132                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1133                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1134                          int_local[i]+=(f_0+f_1)*w_0;                          int_local[i]+=(f0+f1)*w0;
1135                      }  /* end of component loop i */                      }  // end of component loop i
1136                  } /* end of k0 loop */                  } // end of k0 loop
1137              }              }
1138    
1139              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1140  #pragma omp for nowait  #pragma omp for nowait
1141                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1142                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1143                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1144                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1145                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1146                          int_local[i]+=(f_0+f_1)*w_0;                          int_local[i]+=(f0+f1)*w0;
1147                      }  /* end of component loop i */                      }  // end of component loop i
1148                  } /* end of k0 loop */                  } // end of k0 loop
1149              }              }
   
1150  #pragma omp critical  #pragma omp critical
1151              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1152                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1153          } // end of parallel section          } // end of parallel section
1154      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
1155        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1156  #pragma omp parallel  #pragma omp parallel
1157          {          {
1158              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1159              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1160  #pragma omp for nowait  #pragma omp for nowait
1161                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1162                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1163                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1164                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*m_dx[1];
1165                      }  /* end of component loop i */                      }
1166                  } /* end of k1 loop */                  }
1167              }              }
1168    
1169              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1170  #pragma omp for nowait  #pragma omp for nowait
1171                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1172                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1173                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1174                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*m_dx[1];
1175                      }  /* end of component loop i */                      }
1176                  } /* end of k1 loop */                  }
1177              }              }
1178    
1179              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1180  #pragma omp for nowait  #pragma omp for nowait
1181                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1182                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1183                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1184                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*m_dx[0];
1185                      }  /* end of component loop i */                      }
1186                  } /* end of k0 loop */                  }
1187              }              }
1188    
1189              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1190  #pragma omp for nowait  #pragma omp for nowait
1191                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1192                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1193                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1194                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*m_dx[0];
1195                      }  /* end of component loop i */                      }
1196                  } /* end of k0 loop */                  }
1197              }              }
1198    
1199  #pragma omp critical  #pragma omp critical
1200              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1201                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1202          } // end of parallel section          } // end of parallel section
1203      } else {      } // function space selector
         stringstream msg;  
         msg << "setToIntegrals() not implemented for "  
             << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
1204  }  }
1205    
1206  void Rectangle::setToNormal(escript::Data& out) const  //protected
1207    dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
1208  {  {
1209      if (out.getFunctionSpace().getTypeCode() == FaceElements) {      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1210  #pragma omp parallel      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1211          {      const int x=node%nDOF0;
1212              if (m_faceOffset[0] > -1) {      const int y=node/nDOF0;
1213  #pragma omp for nowait      dim_t num=0;
1214                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {      // loop through potential neighbours and add to index if positions are
1215                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);      // within bounds
1216                      // set vector at two quadrature points      for (int i1=-1; i1<2; i1++) {
1217                      *o++ = -1.;          for (int i0=-1; i0<2; i0++) {
1218                      *o++ = 0.;              // skip node itself
1219                      *o++ = -1.;              if (i0==0 && i1==0)
1220                      *o = 0.;                  continue;
1221                  }              // location of neighbour node
1222              }              const int nx=x+i0;
1223                const int ny=y+i1;
1224              if (m_faceOffset[1] > -1) {              if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {
1225  #pragma omp for nowait                  index.push_back(ny*nDOF0+nx);
1226                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  num++;
                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);  
                     // set vector at two quadrature points  
                     *o++ = 1.;  
                     *o++ = 0.;  
                     *o++ = 1.;  
                     *o = 0.;  
                 }  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);  
                     // set vector at two quadrature points  
                     *o++ = 0.;  
                     *o++ = -1.;  
                     *o++ = 0.;  
                     *o = -1.;  
                 }  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);  
                     // set vector at two quadrature points  
                     *o++ = 0.;  
                     *o++ = 1.;  
                     *o++ = 0.;  
                     *o = 1.;  
                 }  
             }  
         } // end of parallel section  
     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
 #pragma omp parallel  
         {  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);  
                     *o++ = -1.;  
                     *o = 0.;  
                 }  
1227              }              }
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);  
                     *o++ = 1.;  
                     *o = 0.;  
                 }  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);  
                     *o++ = 0.;  
                     *o = -1.;  
                 }  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);  
                     *o++ = 0.;  
                     *o = 1.;  
                 }  
             }  
         } // end of parallel section  
   
     } else {  
         stringstream msg;  
         msg << "setToNormal() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
 }  
   
 void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,  
         escript::Data& rhs, const escript::Data& A, const escript::Data& B,  
         const escript::Data& C, const escript::Data& D,  
         const escript::Data& X, const escript::Data& Y,  
         const escript::Data& d, const escript::Data& y,  
         const escript::Data& d_contact, const escript::Data& y_contact,  
         const escript::Data& d_dirac, const escript::Data& y_dirac) const  
 {  
     throw RipleyException("addPDEToSystem() not implemented");  
 }  
   
 Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  
                                                 bool reducedColOrder) const  
 {  
     if (reducedRowOrder || reducedColOrder)  
         throw RipleyException("getPattern() not implemented for reduced order");  
   
     // connector  
     RankVector neighbour;  
     IndexVector offsetInShared(1,0);  
     IndexVector sendShared, recvShared;  
     const IndexVector faces=getNumFacesPerBoundary();  
     const index_t left = (m_offset0==0 ? 0 : 1);  
     const index_t bottom = (m_offset1==0 ? 0 : 1);  
     // corner node from bottom-left  
     if (faces[0] == 0 && faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX-1);  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0+left]);  
         recvShared.push_back(m_nodeId[0]);  
     }  
     // bottom edge  
     if (faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX);  
         offsetInShared.push_back(offsetInShared.back()+m_N0-left);  
         for (dim_t i=left; i<m_N0; i++) {  
             // easy case, we know the neighbour id's  
             sendShared.push_back(m_nodeId[m_N0+i]);  
             recvShared.push_back(m_nodeId[i]);  
         }  
     }  
     // corner node from bottom-right  
     if (faces[1] == 0 && faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX+1);  
         const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour.back()/m_NX == 0 ? m_N1 : m_N1-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[(bottom+1)*m_N0-1]);  
         recvShared.push_back(first+N0*(N1-1));  
     }  
     // left edge  
     if (faces[0] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-1);  
         offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);  
         for (dim_t i=bottom; i<m_N1; i++) {  
             // easy case, we know the neighbour id's  
             sendShared.push_back(m_nodeId[i*m_N0+1]);  
             recvShared.push_back(m_nodeId[i*m_N0]);  
         }  
     }  
     // right edge  
     if (faces[1] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+1);  
         const index_t rightN0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);  
         for (dim_t i=bottom; i<m_N1; i++) {  
             sendShared.push_back(m_nodeId[(i+1)*m_N0-1]);  
             recvShared.push_back(first+rightN0*(i-bottom));  
         }  
     }  
     // corner node from top-left  
     if (faces[0] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX-1);  
         const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+left]);  
         recvShared.push_back(first+N0-1);  
     }  
     // top edge  
     if (faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+m_N0-left);  
         for (dim_t i=left; i<m_N0; i++) {  
             sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+i]);  
             recvShared.push_back(first+i-left);  
         }  
     }  
     // corner node from top-right  
     if (faces[1] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX+1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0*m_N1-1]);  
         recvShared.push_back(first);  
     }  
     const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];  
     /*  
     cout << "--- rcv_shcomp ---" << endl;  
     cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;  
     for (size_t i=0; i<neighbour.size(); i++) {  
         cout << "neighbor[" << i << "]=" << neighbour[i]  
             << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;  
     }  
     for (size_t i=0; i<recvShared.size(); i++) {  
         cout << "shared[" << i << "]=" << recvShared[i] << endl;  
     }  
     cout << "--- snd_shcomp ---" << endl;  
     for (size_t i=0; i<sendShared.size(); i++) {  
         cout << "shared[" << i << "]=" << sendShared[i] << endl;  
     }  
     */  
   
     Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &sendShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &recvShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);  
     Paso_SharedComponents_free(snd_shcomp);  
     Paso_SharedComponents_free(rcv_shcomp);  
   
     // create patterns  
     dim_t M, N;  
     IndexVector ptr(1,0);  
     IndexVector index;  
   
     // main pattern  
     for (index_t i=0; i<numDOF; i++) {  
         // always add the node itself  
         index.push_back(i);  
         int num=insertNeighbours(index, i);  
         ptr.push_back(ptr.back()+num+1);  
     }  
     M=N=ptr.size()-1;  
     // paso will manage the memory  
     index_t* indexC = MEMALLOC(index.size(),index_t);  
     index_t* ptrC = MEMALLOC(ptr.size(), index_t);  
     copy(index.begin(), index.end(), indexC);  
     copy(ptr.begin(), ptr.end(), ptrC);  
     Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
             M, N, ptrC, indexC);  
   
     /*  
     cout << "--- main_pattern ---" << endl;  
     cout << "M=" << M << ", N=" << N << endl;  
     for (size_t i=0; i<ptr.size(); i++) {  
         cout << "ptr[" << i << "]=" << ptr[i] << endl;  
     }  
     for (size_t i=0; i<index.size(); i++) {  
         cout << "index[" << i << "]=" << index[i] << endl;  
     }  
     */  
   
     ptr.clear();  
     index.clear();  
   
     // column & row couple patterns  
     Paso_Pattern *colCouplePattern, *rowCouplePattern;  
     generateCouplePatterns(&colCouplePattern, &rowCouplePattern);  
   
     // allocate paso distribution  
     Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,  
             const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);  
   
     Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(  
             MATRIX_FORMAT_DEFAULT, distribution, distribution,  
             mainPattern, colCouplePattern, rowCouplePattern,  
             connector, connector);  
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colCouplePattern);  
     Paso_Pattern_free(rowCouplePattern);  
     Paso_Distribution_free(distribution);  
     return pattern;  
 }  
   
 void Rectangle::Print_Mesh_Info(const bool full) const  
 {  
     RipleyDomain::Print_Mesh_Info(full);  
     if (full) {  
         cout << "     Id  Coordinates" << endl;  
         cout.precision(15);  
         cout.setf(ios::scientific, ios::floatfield);  
         pair<double,double> xdx = getFirstCoordAndSpacing(0);  
         pair<double,double> ydy = getFirstCoordAndSpacing(1);  
         for (index_t i=0; i < getNumNodes(); i++) {  
             cout << "  " << setw(5) << m_nodeId[i]  
                 << "  " << xdx.first+(i%m_N0)*xdx.second  
                 << "  " << ydy.first+(i/m_N0)*ydy.second << endl;  
1228          }          }
1229      }      }
 }  
1230    
1231  IndexVector Rectangle::getNumNodesPerDim() const      return num;
 {  
     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;  
 }  
   
 pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  
 {  
     if (dim==0) {  
         return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);  
     } else if (dim==1) {  
         return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);  
     }  
     throw RipleyException("getFirstCoordAndSpacing(): invalid argument");  
1232  }  }
1233    
1234  //protected  //protected
1235  dim_t Rectangle::getNumDOF() const  void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const
1236  {  {
1237      return m_nodeDistribution[m_mpiInfo->rank+1]      const dim_t numComp = in.getDataPointSize();
1238          -m_nodeDistribution[m_mpiInfo->rank];      out.requireWrite();
 }  
1239    
1240  //protected      const index_t left = (m_offset[0]==0 ? 0 : 1);
1241  dim_t Rectangle::getNumFaceElements() const      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1242  {      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1243      const IndexVector faces = getNumFacesPerBoundary();      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1244      dim_t n=0;  #pragma omp parallel for
1245      for (size_t i=0; i<faces.size(); i++)      for (index_t i=0; i<nDOF1; i++) {
1246          n+=faces[i];          for (index_t j=0; j<nDOF0; j++) {
1247      return n;              const index_t n=j+left+(i+bottom)*m_NN[0];
1248                const double* src=in.getSampleDataRO(n);
1249                copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
1250            }
1251        }
1252  }  }
1253    
1254  //protected  //protected
1255  void Rectangle::assembleCoordinates(escript::Data& arg) const  void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const
1256  {  {
1257      escriptDataC x = arg.getDataC();      const dim_t numComp = in.getDataPointSize();
1258      int numDim = m_numDim;      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
1259      if (!isDataPointShapeEqual(&x, 1, &numDim))      in.requireWrite();
1260          throw RipleyException("setToX: Invalid Data object shape");      Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));
1261      if (!numSamplesEqual(&x, 1, getNumNodes()))  
1262          throw RipleyException("setToX: Illegal number of samples in Data object");      const dim_t numDOF = getNumDOF();
1263        out.requireWrite();
1264        const double* buffer = Paso_Coupler_finishCollect(coupler);
1265    
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
     arg.requireWrite();  
1266  #pragma omp parallel for  #pragma omp parallel for
1267      for (dim_t i1 = 0; i1 < m_N1; i1++) {      for (index_t i=0; i<getNumNodes(); i++) {
1268          for (dim_t i0 = 0; i0 < m_N0; i0++) {          const double* src=(m_dofMap[i]<numDOF ?
1269              double* point = arg.getSampleDataRW(i0+m_N0*i1);                  in.getSampleDataRO(m_dofMap[i])
1270              point[0] = xdx.first+i0*xdx.second;                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1271              point[1] = ydy.first+i1*ydy.second;          copy(src, src+numComp, out.getSampleDataRW(i));
         }  
1272      }      }
1273        Paso_Coupler_free(coupler);
1274  }  }
1275    
1276  //private  //private
1277  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
1278  {  {
1279      // identifiers are ordered from left to right, bottom to top on each rank,      // degrees of freedom are numbered from left to right, bottom to top in
1280      // except for the shared nodes which are owned by the rank below / to the      // each rank, continuing on the next rank (ranks also go left-right,
1281      // left of the current rank      // bottom-top).
1282        // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
1283        // helps when writing out data rank after rank.
1284    
1285      // build node distribution vector first.      // build node distribution vector first.
1286      // m_nodeDistribution[i] is the first node id on rank i, that is      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
1287      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // constant for all ranks in this implementation
1288      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1289      m_nodeDistribution[1]=getNumNodes();      const dim_t numDOF=getNumDOF();
1290      for (dim_t k=1; k<m_mpiInfo->size-1; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
1291          const index_t x=k%m_NX;          m_nodeDistribution[k]=k*numDOF;
         const index_t y=k/m_NX;  
         index_t numNodes=getNumNodes();  
         if (x>0)  
             numNodes-=m_N1;  
         if (y>0)  
             numNodes-=m_N0;  
         if (x>0 && y>0)  
             numNodes++; // subtracted corner twice -> fix that  
         m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;  
1292      }      }
1293      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
   
     m_dofId.resize(getNumDOF());  
1294      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
1295        m_dofId.resize(numDOF);
1296        m_elementId.resize(getNumElements());
1297    
1298      // the bottom row and left column are not owned by this rank so the      // populate face element counts
1299      // identifiers need to be computed accordingly      //left
1300      const index_t left = (m_offset0==0 ? 0 : 1);      if (m_offset[0]==0)
1301      const index_t bottom = (m_offset1==0 ? 0 : 1);          m_faceCount[0]=m_NE[1];
1302      if (left>0) {      else
1303          const int neighbour=m_mpiInfo->rank-1;          m_faceCount[0]=0;
1304          const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      //right
1305  #pragma omp parallel for      if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
1306          for (dim_t i1=bottom; i1<m_N1; i1++) {          m_faceCount[1]=m_NE[1];
1307              m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]      else
1308                  + (i1-bottom+1)*leftN0-1;          m_faceCount[1]=0;
1309        //bottom
1310        if (m_offset[1]==0)
1311            m_faceCount[2]=m_NE[0];
1312        else
1313            m_faceCount[2]=0;
1314        //top
1315        if (m_mpiInfo->rank/m_NX[0]==m_NX[1]-1)
1316            m_faceCount[3]=m_NE[0];
1317        else
1318            m_faceCount[3]=0;
1319    
1320        m_faceId.resize(getNumFaceElements());
1321    
1322        const index_t left = (m_offset[0]==0 ? 0 : 1);
1323        const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1324        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1325        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1326    
1327    #define globalNodeId(x,y) \
1328        ((m_offset[0]+x)/nDOF0)*nDOF0*nDOF1+(m_offset[0]+x)%nDOF0 \
1329        + ((m_offset[1]+y)/nDOF1)*nDOF0*nDOF1*m_NX[0]+((m_offset[1]+y)%nDOF1)*nDOF0
1330    
1331        // set corner id's outside the parallel region
1332        m_nodeId[0] = globalNodeId(0, 0);
1333        m_nodeId[m_NN[0]-1] = globalNodeId(m_NN[0]-1, 0);
1334        m_nodeId[m_NN[0]*(m_NN[1]-1)] = globalNodeId(0, m_NN[1]-1);
1335        m_nodeId[m_NN[0]*m_NN[1]-1] = globalNodeId(m_NN[0]-1,m_NN[1]-1);
1336    #undef globalNodeId
1337    
1338    #pragma omp parallel
1339        {
1340            // populate degrees of freedom and own nodes (identical id)
1341    #pragma omp for nowait
1342            for (dim_t i=0; i<nDOF1; i++) {
1343                for (dim_t j=0; j<nDOF0; j++) {
1344                    const index_t nodeIdx=j+left+(i+bottom)*m_NN[0];
1345                    const index_t dofIdx=j+i*nDOF0;
1346                    m_dofId[dofIdx] = m_nodeId[nodeIdx]
1347                        = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
1348                }
1349          }          }
1350      }  
1351      if (bottom>0) {          // populate the rest of the nodes (shared with other ranks)
1352          const int neighbour=m_mpiInfo->rank-m_NX;          if (m_faceCount[0]==0) { // left column
1353          const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  #pragma omp for nowait
1354          const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);              for (dim_t i=0; i<nDOF1; i++) {
1355  #pragma omp parallel for                  const index_t nodeIdx=(i+bottom)*m_NN[0];
1356          for (dim_t i0=left; i0<m_N0; i0++) {                  const index_t dofId=(i+1)*nDOF0-1;
1357              m_nodeId[i0]=m_nodeDistribution[neighbour]                  m_nodeId[nodeIdx]
1358                  + (bottomN1-1)*bottomN0 + i0 - left;                      = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
1359                }
1360            }
1361            if (m_faceCount[1]==0) { // right column
1362    #pragma omp for nowait
1363                for (dim_t i=0; i<nDOF1; i++) {
1364                    const index_t nodeIdx=(i+bottom+1)*m_NN[0]-1;
1365                    const index_t dofId=i*nDOF0;
1366                    m_nodeId[nodeIdx]
1367                        = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
1368                }
1369            }
1370            if (m_faceCount[2]==0) { // bottom row
1371    #pragma omp for nowait
1372                for (dim_t i=0; i<nDOF0; i++) {
1373                    const index_t nodeIdx=i+left;
1374                    const index_t dofId=nDOF0*(nDOF1-1)+i;
1375                    m_nodeId[nodeIdx]
1376                        = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
1377                }
1378            }
1379            if (m_faceCount[3]==0) { // top row
1380    #pragma omp for nowait
1381                for (dim_t i=0; i<nDOF0; i++) {
1382                    const index_t nodeIdx=m_NN[0]*(m_NN[1]-1)+i+left;
1383                    const index_t dofId=i;
1384                    m_nodeId[nodeIdx]
1385                        = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
1386                }
1387          }          }
     }  
     if (left>0 && bottom>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX-1;  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);  
         m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1-1;  
     }  
1388    
1389      // the rest of the id's are contiguous          // populate element id's
1390      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];  #pragma omp for nowait
1391  #pragma omp parallel for          for (dim_t i1=0; i1<m_NE[1]; i1++) {
1392      for (dim_t i1=bottom; i1<m_N1; i1++) {              for (dim_t i0=0; i0<m_NE[0]; i0++) {
1393          for (dim_t i0=left; i0<m_N0; i0++) {                  m_elementId[i0+i1*m_NE[0]]=(m_offset[1]+i1)*m_gNE[0]+m_offset[0]+i0;
1394              const index_t idx=i0-left+(i1-bottom)*(m_N0-left);              }
             m_nodeId[i0+i1*m_N0] = firstId+idx;  
             m_dofId[idx] = firstId+idx;  
1395          }          }
1396      }  
1397            // face elements
1398    #pragma omp for
1399            for (dim_t k=0; k<getNumFaceElements(); k++)
1400                m_faceId[k]=k;
1401        } // end parallel section
1402    
1403      m_nodeTags.assign(getNumNodes(), 0);      m_nodeTags.assign(getNumNodes(), 0);
1404      updateTagsInUse(Nodes);      updateTagsInUse(Nodes);
1405    
     // elements  
     m_elementId.resize(getNumElements());  
 #pragma omp parallel for  
     for (dim_t k=0; k<getNumElements(); k++) {  
         m_elementId[k]=k;  
     }  
1406      m_elementTags.assign(getNumElements(), 0);      m_elementTags.assign(getNumElements(), 0);
1407      updateTagsInUse(Elements);      updateTagsInUse(Elements);
1408    
     // face element id's  
     m_faceId.resize(getNumFaceElements());  
 #pragma omp parallel for  
     for (dim_t k=0; k<getNumFaceElements(); k++) {  
         m_faceId[k]=k;  
     }  
   
1409      // generate face offset vector and set face tags      // generate face offset vector and set face tags
     const IndexVector facesPerEdge = getNumFacesPerBoundary();  
1410      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
1411      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };      const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
1412      m_faceOffset.assign(facesPerEdge.size(), -1);      m_faceOffset.assign(4, -1);
1413      m_faceTags.clear();      m_faceTags.clear();
1414      index_t offset=0;      index_t offset=0;
1415      for (size_t i=0; i<facesPerEdge.size(); i++) {      for (size_t i=0; i<4; i++) {
1416          if (facesPerEdge[i]>0) {          if (m_faceCount[i]>0) {
1417              m_faceOffset[i]=offset;              m_faceOffset[i]=offset;
1418              offset+=facesPerEdge[i];              offset+=m_faceCount[i];
1419              m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);              m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
1420          }          }
1421      }      }
1422      setTagMap("left", LEFT);      setTagMap("left", LEFT);
# Line 1121  void Rectangle::populateSampleIds() Line 1427  void Rectangle::populateSampleIds()
1427  }  }
1428    
1429  //private  //private
1430  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const  void Rectangle::createPattern()
1431  {  {
1432      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1433      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1434      const int x=node%myN0;      const index_t left = (m_offset[0]==0 ? 0 : 1);
1435      const int y=node/myN0;      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1436      int num=0;  
1437      if (y>0) {      // populate node->DOF mapping with own degrees of freedom.
1438          if (x>0) {      // The rest is assigned in the loop further down
1439              // bottom-left      m_dofMap.assign(getNumNodes(), 0);
1440              index.push_back(node-myN0-1);  #pragma omp parallel for
1441              num++;      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1442          }          for (index_t j=left; j<left+nDOF0; j++) {
1443          // bottom              m_dofMap[i*m_NN[0]+j]=(i-bottom)*nDOF0+j-left;
1444          index.push_back(node-myN0);          }
         num++;  
         if (x<myN0-1) {  
             // bottom-right  
             index.push_back(node-myN0+1);  
             num++;  
         }  
     }  
     if (x<myN0-1) {  
         // right  
         index.push_back(node+1);  
         num++;  
         if (y<myN1-1) {  
             // top-right  
             index.push_back(node+myN0+1);  
             num++;  
         }  
     }  
     if (y<myN1-1) {  
         // top  
         index.push_back(node+myN0);  
         num++;  
         if (x>0) {  
             // top-left  
             index.push_back(node+myN0-1);  
             num++;  
         }  
     }  
     if (x>0) {  
         // left  
         index.push_back(node-1);  
         num++;  
1445      }      }
1446    
1447      return num;      // build list of shared components and neighbours by looping through
1448  }      // all potential neighbouring ranks and checking if positions are
1449        // within bounds
1450  //private      const dim_t numDOF=nDOF0*nDOF1;
1451  void Rectangle::generateCouplePatterns(Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const      vector<IndexVector> colIndices(numDOF); // for the couple blocks
1452  {      RankVector neighbour;
1453      IndexVector ptr(1,0);      IndexVector offsetInShared(1,0);
1454      IndexVector index;      IndexVector sendShared, recvShared;
1455      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);      int numShared=0;
1456      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);      const int x=m_mpiInfo->rank%m_NX[0];
1457      const IndexVector faces=getNumFacesPerBoundary();      const int y=m_mpiInfo->rank/m_NX[0];
1458        for (int i1=-1; i1<2; i1++) {
1459      // bottom edge          for (int i0=-1; i0<2; i0++) {
1460      dim_t n=0;              // skip this rank
1461      if (faces[0] == 0) {              if (i0==0 && i1==0)
1462          index.push_back(2*(myN0+myN1+1));                  continue;
1463          index.push_back(2*(myN0+myN1+1)+1);              // location of neighbour rank
1464          n+=2;              const int nx=x+i0;
1465          if (faces[2] == 0) {              const int ny=y+i1;
1466              index.push_back(0);              if (nx>=0 && ny>=0 && nx<m_NX[0] && ny<m_NX[1]) {
1467              index.push_back(1);                  neighbour.push_back(ny*m_NX[0]+nx);
1468              index.push_back(2);                  if (i0==0) {
1469              n+=3;                      // sharing top or bottom edge
1470          }                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1471      } else if (faces[2] == 0) {                      const int firstNode=(i1==-1 ? left : m_NN[0]*(m_NN[1]-1)+left);
1472          index.push_back(1);                      offsetInShared.push_back(offsetInShared.back()+nDOF0);
1473          index.push_back(2);                      for (dim_t i=0; i<nDOF0; i++, numShared++) {
1474          n+=2;                          sendShared.push_back(firstDOF+i);
1475      }                          recvShared.push_back(numDOF+numShared);
1476      // n=neighbours of bottom-left corner node                          if (i>0)
1477      ptr.push_back(ptr.back()+n);                              colIndices[firstDOF+i-1].push_back(numShared);
1478      n=0;                          colIndices[firstDOF+i].push_back(numShared);
1479      if (faces[2] == 0) {                          if (i<nDOF0-1)
1480          for (dim_t i=1; i<myN0-1; i++) {                              colIndices[firstDOF+i+1].push_back(numShared);
1481              index.push_back(i);                          m_dofMap[firstNode+i]=numDOF+numShared;
1482              index.push_back(i+1);                      }
1483              index.push_back(i+2);                  } else if (i1==0) {
1484              ptr.push_back(ptr.back()+3);                      // sharing left or right edge
1485          }                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1486          index.push_back(myN0-1);                      const int firstNode=(i0==-1 ? bottom*m_NN[0] : (bottom+1)*m_NN[0]-1);
1487          index.push_back(myN0);                      offsetInShared.push_back(offsetInShared.back()+nDOF1);
1488          n+=2;                      for (dim_t i=0; i<nDOF1; i++, numShared++) {
1489          if (faces[1] == 0) {                          sendShared.push_back(firstDOF+i*nDOF0);
1490              index.push_back(myN0+1);                          recvShared.push_back(numDOF+numShared);
1491              index.push_back(myN0+2);                          if (i>0)
1492              index.push_back(myN0+3);                              colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);
1493              n+=3;                          colIndices[firstDOF+i*nDOF0].push_back(numShared);
1494          }                          if (i<nDOF1-1)
1495      } else {                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1496          for (dim_t i=1; i<myN0-1; i++) {                          m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
1497              ptr.push_back(ptr.back());                      }
1498          }                  } else {
1499          if (faces[1] == 0) {                      // sharing a node
1500              index.push_back(myN0+2);                      const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1501              index.push_back(myN0+3);                      const int node=(i0+1)/2*(m_NN[0]-1)+(i1+1)/2*m_NN[0]*(m_NN[1]-1);
1502              n+=2;                      offsetInShared.push_back(offsetInShared.back()+1);
1503          }                      sendShared.push_back(dof);
1504      }                      recvShared.push_back(numDOF+numShared);
1505      // n=neighbours of bottom-right corner node                      colIndices[dof].push_back(numShared);
1506      ptr.push_back(ptr.back()+n);                      m_dofMap[node]=numDOF+numShared;
1507                        ++numShared;
1508      // 2nd row to 2nd last row                  }
1509      for (dim_t i=1; i<myN1-1; i++) {              }
         // left edge  
         if (faces[0] == 0) {  
             index.push_back(2*(myN0+myN1+2)-i);  
             index.push_back(2*(myN0+myN1+2)-i-1);  
             index.push_back(2*(myN0+myN1+2)-i-2);  
             ptr.push_back(ptr.back()+3);  
         } else {  
             ptr.push_back(ptr.back());  
         }  
         for (dim_t j=1; j<myN0-1; j++) {  
             ptr.push_back(ptr.back());  
         }  
         // right edge  
         if (faces[1] == 0) {  
             index.push_back(myN0+i+1);  
             index.push_back(myN0+i+2);  
             index.push_back(myN0+i+3);  
             ptr.push_back(ptr.back()+3);  
         } else {  
             ptr.push_back(ptr.back());  
         }  
     }  
   
     // top edge  
     n=0;  
     if (faces[0] == 0) {  
         index.push_back(2*myN0+myN1+5);  
         index.push_back(2*myN0+myN1+4);  
         n+=2;  
         if (faces[3] == 0) {  
             index.push_back(2*myN0+myN1+3);  
             index.push_back(2*myN0+myN1+2);  
             index.push_back(2*myN0+myN1+1);  
             n+=3;  
         }  
     } else if (faces[3] == 0) {  
         index.push_back(2*myN0+myN1+2);  
         index.push_back(2*myN0+myN1+1);  
         n+=2;  
     }  
     // n=neighbours of top-left corner node  
     ptr.push_back(ptr.back()+n);  
     n=0;  
     if (faces[3] == 0) {  
         for (dim_t i=1; i<myN0-1; i++) {  
             index.push_back(2*myN0+myN1+i+1);  
             index.push_back(2*myN0+myN1+i);  
             index.push_back(2*myN0+myN1+i-1);  
             ptr.push_back(ptr.back()+3);  
         }  
         index.push_back(myN0+myN1+4);  
         index.push_back(myN0+myN1+3);  
         n+=2;  
         if (faces[1] == 0) {  
             index.push_back(myN0+myN1+2);  
             index.push_back(myN0+myN1+1);  
             index.push_back(myN0+myN1);  
             n+=3;  
         }  
     } else {  
         for (dim_t i=1; i<myN0-1; i++) {  
             ptr.push_back(ptr.back());  
         }  
         if (faces[1] == 0) {  
             index.push_back(myN0+myN1+1);  
             index.push_back(myN0+myN1);  
             n+=2;  
         }  
     }  
     // n=neighbours of top-right corner node  
     ptr.push_back(ptr.back()+n);  
   
     dim_t M=ptr.size()-1;  
     map<index_t,index_t> idMap;  
     dim_t N=0;  
     for (IndexVector::iterator it=index.begin(); it!=index.end(); it++) {  
         if (idMap.find(*it)==idMap.end()) {  
             idMap[*it]=N;  
             *it=N++;  
         } else {  
             *it=idMap[*it];  
1510          }          }
1511      }      }
1512    
1513        // create connector
1514        Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1515                numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
1516                &offsetInShared[0], 1, 0, m_mpiInfo);
1517        Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
1518                numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
1519                &offsetInShared[0], 1, 0, m_mpiInfo);
1520        m_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1521        Paso_SharedComponents_free(snd_shcomp);
1522        Paso_SharedComponents_free(rcv_shcomp);
1523    
1524        // create main and couple blocks
1525        Paso_Pattern *mainPattern = createMainPattern();
1526        Paso_Pattern *colPattern, *rowPattern;
1527        createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
1528    
1529        // allocate paso distribution
1530        Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
1531                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
1532    
1533        // finally create the system matrix
1534        m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
1535                distribution, distribution, mainPattern, colPattern, rowPattern,
1536                m_connector, m_connector);
1537    
1538        Paso_Distribution_free(distribution);
1539    
1540        // useful debug output
1541      /*      /*
1542      cout << "--- colCouple_pattern ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
1543      cout << "M=" << M << ", N=" << N << endl;      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
1544      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<neighbour.size(); i++) {
1545          cout << "ptr[" << i << "]=" << ptr[i] << endl;          cout << "neighbor[" << i << "]=" << neighbour[i]
1546                << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
1547      }      }
1548      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<recvShared.size(); i++) {
1549          cout << "index[" << i << "]=" << index[i] << endl;          cout << "shared[" << i << "]=" << recvShared[i] << endl;
1550        }
1551        cout << "--- snd_shcomp ---" << endl;
1552        for (size_t i=0; i<sendShared.size(); i++) {
1553            cout << "shared[" << i << "]=" << sendShared[i] << endl;
1554        }
1555        cout << "--- dofMap ---" << endl;
1556        for (size_t i=0; i<m_dofMap.size(); i++) {
1557            cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
1558        }
1559        cout << "--- colIndices ---" << endl;
1560        for (size_t i=0; i<colIndices.size(); i++) {
1561            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
1562      }      }
1563      */      */
1564    
1565      // now build the row couple pattern      /*
1566      IndexVector ptr2(1,0);      cout << "--- main_pattern ---" << endl;
1567      IndexVector index2;      cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
1568      for (dim_t id=0; id<N; id++) {      for (size_t i=0; i<mainPattern->numOutput+1; i++) {
1569          n=0;          cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
1570          for (dim_t i=0; i<M; i++) {      }
1571              for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {      for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
1572                  if (index[j]==id) {          cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
                     index2.push_back(i);  
                     n++;  
                     break;  
                 }  
             }  
         }  
         ptr2.push_back(ptr2.back()+n);  
1573      }      }
1574        */
1575    
1576        /*
1577        cout << "--- colCouple_pattern ---" << endl;
1578        cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
1579        for (size_t i=0; i<colPattern->numOutput+1; i++) {
1580            cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
1581        }
1582        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
1583            cout << "index[" << i << "]=" << colPattern->index[i] << endl;
1584        }
1585        */
1586    
1587      /*      /*
1588      cout << "--- rowCouple_pattern ---" << endl;      cout << "--- rowCouple_pattern ---" << endl;
1589      cout << "M=" << N << ", N=" << M << endl;      cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
1590      for (size_t i=0; i<ptr2.size(); i++) {      for (size_t i=0; i<rowPattern->numOutput+1; i++) {
1591          cout << "ptr[" << i << "]=" << ptr2[i] << endl;          cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
1592      }      }
1593      for (size_t i=0; i<index2.size(); i++) {      for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
1594          cout << "index[" << i << "]=" << index2[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
1595      }      }
1596      */      */
1597    
1598      // paso will manage the memory      Paso_Pattern_free(mainPattern);
1599      index_t* indexC = MEMALLOC(index.size(), index_t);      Paso_Pattern_free(colPattern);
1600      index_t* ptrC = MEMALLOC(ptr.size(), index_t);      Paso_Pattern_free(rowPattern);
1601      copy(index.begin(), index.end(), indexC);  }
1602      copy(ptr.begin(), ptr.end(), ptrC);  
1603      *colPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);  //private
1604    void Rectangle::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
1605      // paso will manage the memory           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
1606      indexC = MEMALLOC(index2.size(), index_t);           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
1607      ptrC = MEMALLOC(ptr2.size(), index_t);  {
1608      copy(index2.begin(), index2.end(), indexC);      IndexVector rowIndex;
1609      copy(ptr2.begin(), ptr2.end(), ptrC);      rowIndex.push_back(m_dofMap[firstNode]);
1610      *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);      rowIndex.push_back(m_dofMap[firstNode+1]);
1611        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]]);
1612        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]+1]);
1613        if (addF) {
1614            double *F_p=F.getSampleDataRW(0);
1615            for (index_t i=0; i<rowIndex.size(); i++) {
1616                if (rowIndex[i]<getNumDOF()) {
1617                    for (index_t eq=0; eq<nEq; eq++) {
1618                        F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
1619                    }
1620                }
1621            }
1622        }
1623        if (addS) {
1624            addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
1625        }
1626  }  }
1627    
1628  //protected  //protected
# Line 1376  void Rectangle::interpolateNodesOnElemen Line 1631  void Rectangle::interpolateNodesOnElemen
1631  {  {
1632      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1633      if (reduced) {      if (reduced) {
1634          /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */          out.requireWrite();
1635          const double c0 = .25;          const double c0 = 0.25;
1636  #pragma omp parallel for  #pragma omp parallel
1637          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1638              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1639                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_01(numComp);
1640                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1641                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              vector<double> f_11(numComp);
1642                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1643                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1644                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1645                      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));
1646                  } /* end of component loop i */                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_NN[0])), numComp*sizeof(double));
1647              } /* end of k0 loop */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_NN[0])), numComp*sizeof(double));
1648          } /* end of k1 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_NN[0])), numComp*sizeof(double));
1649          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1650                        for (index_t i=0; i < numComp; ++i) {
1651                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1652                        } /* end of component loop i */
1653                    } /* end of k0 loop */
1654                } /* end of k1 loop */
1655            } /* end of parallel section */
1656      } else {      } else {
1657          /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */          out.requireWrite();
1658          const double c0 = .16666666666666666667;          const double c0 = 0.16666666666666666667;
1659          const double c1 = .044658198738520451079;          const double c1 = 0.044658198738520451079;
1660          const double c2 = .62200846792814621559;          const double c2 = 0.62200846792814621559;
1661  #pragma omp parallel for  #pragma omp parallel
1662          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1663              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1664                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
1665                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1666                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
1667                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1668                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1669                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1670                      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));
1671                      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));
1672                      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));
1673                      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));
1674                  } /* end of component loop i */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE[0]));
1675              } /* end of k0 loop */                      for (index_t i=0; i < numComp; ++i) {
1676          } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1677          /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */                          o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1678                            o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1679                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1680                        } /* end of component loop i */
1681                    } /* end of k0 loop */
1682                } /* end of k1 loop */
1683            } /* end of parallel section */
1684      }      }
1685  }  }
1686    
# Line 1423  void Rectangle::interpolateNodesOnFaces( Line 1690  void Rectangle::interpolateNodesOnFaces(
1690  {  {
1691      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1692      if (reduced) {      if (reduced) {
1693          const double c0 = .5;          out.requireWrite();
1694  #pragma omp parallel  #pragma omp parallel
1695          {          {
1696              /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */              vector<double> f_00(numComp);
1697                vector<double> f_01(numComp);
1698                vector<double> f_10(numComp);
1699                vector<double> f_11(numComp);
1700              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1701  #pragma omp for nowait  #pragma omp for nowait
1702                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1703                      const register 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));
1704                      const register 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));
1705                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1706                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1707                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = (f_00[i] + f_01[i])/2;
1708                      } /* end of component loop i */                      } /* end of component loop i */
1709                  } /* end of k1 loop */                  } /* end of k1 loop */
1710              } /* end of face 0 */              } /* end of face 0 */
1711              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1712  #pragma omp for nowait  #pragma omp for nowait
1713                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1714                      const register 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));
1715                      const register 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));
1716                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1717                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1718                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = (f_10[i] + f_11[i])/2;
1719                      } /* end of component loop i */                      } /* end of component loop i */
1720                  } /* end of k1 loop */                  } /* end of k1 loop */
1721              } /* end of face 1 */              } /* end of face 1 */
1722              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1723  #pragma omp for nowait  #pragma omp for nowait
1724                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1725                      const register 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));
1726                      const register 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));
1727                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1728                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1729                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = (f_00[i] + f_10[i])/2;
1730                      } /* end of component loop i */                      } /* end of component loop i */
1731                  } /* end of k0 loop */                  } /* end of k0 loop */
1732              } /* end of face 2 */              } /* end of face 2 */
1733              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1734  #pragma omp for nowait  #pragma omp for nowait
1735                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1736                      const register 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));
1737                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_NN[1]-1, m_NN[0])), numComp*sizeof(double));
1738                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1739                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1740                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = (f_01[i] + f_11[i])/2;
1741                      } /* end of component loop i */                      } /* end of component loop i */
1742                  } /* end of k0 loop */                  } /* end of k0 loop */
1743              } /* end of face 3 */              } /* end of face 3 */
1744              /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */          } /* end of parallel section */
         } // end of parallel section  
1745      } else {      } else {
1746            out.requireWrite();
1747          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1748          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1749  #pragma omp parallel  #pragma omp parallel
1750          {          {
1751              /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */              vector<double> f_00(numComp);
1752                vector<double> f_01(numComp);
1753                vector<double> f_10(numComp);
1754                vector<double> f_11(numComp);
1755              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1756  #pragma omp for nowait  #pragma omp for nowait
1757                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1758                      const register 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));
1759                      const register 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));
1760                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1761                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1762                          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];
1763                          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];
1764                      } /* end of component loop i */                      } /* end of component loop i */
1765                  } /* end of k1 loop */                  } /* end of k1 loop */
1766              } /* end of face 0 */              } /* end of face 0 */
1767              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1768  #pragma omp for nowait  #pragma omp for nowait
1769                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1770                      const register 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));
1771                      const register 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));
1772                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1773                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1774                          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];
1775                          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];
1776                      } /* end of component loop i */                      } /* end of component loop i */
1777                  } /* end of k1 loop */                  } /* end of k1 loop */
1778              } /* end of face 1 */              } /* end of face 1 */
1779              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1780  #pragma omp for nowait  #pragma omp for nowait
1781                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1782                      const register 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));
1783                      const register 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));
1784                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1785                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1786                          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];
1787                          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];
1788                      } /* end of component loop i */                      } /* end of component loop i */
1789                  } /* end of k0 loop */                  } /* end of k0 loop */
1790              } /* end of face 2 */              } /* end of face 2 */
1791              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1792  #pragma omp for nowait  #pragma omp for nowait
1793                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1794                      const register 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));
1795                      const register 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));
1796                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1797                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1798                          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];
1799                          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];
1800                      } /* end of component loop i */                      } /* end of component loop i */
1801                  } /* end of k0 loop */                  } /* end of k0 loop */
1802              } /* end of face 3 */              } /* end of face 3 */
1803              /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */          } /* end of parallel section */
         } // end of parallel section  
1804      }      }
1805  }  }
1806    
1807    //protected
1808    void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1809            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1810            const escript::Data& C, const escript::Data& D,
1811            const escript::Data& X, const escript::Data& Y) const
1812    {
1813        const double SQRT3 = 1.73205080756887719318;
1814        const double w1 = 1.0/24.0;
1815        const double w5 = -SQRT3/24 + 1.0/12;
1816        const double w2 = -SQRT3/24 - 1.0/12;
1817        const double w19 = -m_dx[0]/12;
1818        const double w11 = w19*(SQRT3 + 3)/12;
1819        const double w14 = w19*(-SQRT3 + 3)/12;
1820        const double w16 = w19*(5*SQRT3 + 9)/12;
1821        const double w17 = w19*(-5*SQRT3 + 9)/12;
1822        const double w27 = w19*(-SQRT3 - 3)/2;
1823        const double w28 = w19*(SQRT3 - 3)/2;
1824        const double w18 = -m_dx[1]/12;
1825        const double w12 = w18*(5*SQRT3 + 9)/12;
1826        const double w13 = w18*(-5*SQRT3 + 9)/12;
1827        const double w10 = w18*(SQRT3 + 3)/12;
1828        const double w15 = w18*(-SQRT3 + 3)/12;
1829        const double w25 = w18*(-SQRT3 - 3)/2;
1830        const double w26 = w18*(SQRT3 - 3)/2;
1831        const double w22 = m_dx[0]*m_dx[1]/144;
1832        const double w20 = w22*(SQRT3 + 2);
1833        const double w21 = w22*(-SQRT3 + 2);
1834        const double w23 = w22*(4*SQRT3 + 7);
1835        const double w24 = w22*(-4*SQRT3 + 7);
1836        const double w3 = m_dx[0]/(24*m_dx[1]);
1837        const double w7 = w3*(SQRT3 + 2);
1838        const double w8 = w3*(-SQRT3 + 2);
1839        const double w6 = -m_dx[1]/(24*m_dx[0]);
1840        const double w0 = w6*(SQRT3 + 2);
1841        const double w4 = w6*(-SQRT3 + 2);
1842    
1843        rhs.requireWrite();
1844    #pragma omp parallel
1845        {
1846            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1847    #pragma omp for
1848                for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
1849                    for (index_t k0=0; k0<m_NE[0]; ++k0)  {
1850                        bool add_EM_S=false;
1851                        bool add_EM_F=false;
1852                        vector<double> EM_S(4*4, 0);
1853                        vector<double> EM_F(4, 0);
1854                        const index_t e = k0 + m_NE[0]*k1;
1855                        ///////////////
1856                        // process A //
1857                        ///////////////
1858                        if (!A.isEmpty()) {
1859                            add_EM_S = true;
1860                            const double* A_p = const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1861                            if (A.actsExpanded()) {
1862                                const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1863                                const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1864                                const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1865                                const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1866                                const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1867                                const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1868                                const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1869                                const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1870                                const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1871                                const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1872                                const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1873                                const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1874                                const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1875                                const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1876                                const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1877                                const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1878                                const double tmp0 = w3*(A_11_0 + A_11_1 + A_11_2 + A_11_3);
1879                                const double tmp1 = w1*(A_01_0 + A_01_3 - A_10_1 - A_10_2);
1880                                const double tmp2 = w4*(A_00_2 + A_00_3);
1881                                const double tmp3 = w0*(A_00_0 + A_00_1);
1882                                const double tmp4 = w5*(A_01_2 - A_10_3);
1883                                const double tmp5 = w2*(-A_01_1 + A_10_0);
1884                                const double tmp6 = w5*(A_01_3 + A_10_0);
1885                                const double tmp7 = w3*(-A_11_0 - A_11_1 - A_11_2 - A_11_3);
1886                                const double tmp8 = w6*(A_00_0 + A_00_1 + A_00_2 + A_00_3);
1887                                const double tmp9 = w1*(A_01_1 + A_01_2 + A_10_1 + A_10_2);
1888                                const double tmp10 = w2*(-A_01_0 - A_10_3);
1889                                const double tmp11 = w4*(A_00_0 + A_00_1);
1890                                const double tmp12 = w0*(A_00_2 + A_00_3);
1891                                const double tmp13 = w5*(A_01_1 - A_10_0);
1892                                const double tmp14 = w2*(-A_01_2 + A_10_3);
1893                                const double tmp15 = w7*(A_11_0 + A_11_2);
1894                                const double tmp16 = w4*(-A_00_2 - A_00_3);
1895                                const double tmp17 = w0*(-A_00_0 - A_00_1);
1896                                const double tmp18 = w5*(A_01_3 + A_10_3);
1897                                const double tmp19 = w8*(A_11_1 + A_11_3);
1898                                const double tmp20 = w2*(-A_01_0 - A_10_0);
1899                                const double tmp21 = w7*(A_11_1 + A_11_3);
1900                                const double tmp22 = w4*(-A_00_0 - A_00_1);
1901                                const double tmp23 = w0*(-A_00_2 - A_00_3);
1902                                const double tmp24 = w5*(A_01_0 + A_10_0);
1903                                const double tmp25 = w8*(A_11_0 + A_11_2);
1904                                const double tmp26 = w2*(-A_01_3 - A_10_3);
1905                                const double tmp27 = w5*(-A_01_1 - A_10_2);
1906                                const double tmp28 = w1*(-A_01_0 - A_01_3 - A_10_0 - A_10_3);
1907                                const double tmp29 = w2*(A_01_2 + A_10_1);
1908                                const double tmp30 = w7*(-A_11_1 - A_11_3);
1909                                const double tmp31 = w1*(-A_01_1 - A_01_2 + A_10_0 + A_10_3);
1910                                const double tmp32 = w5*(-A_01_0 + A_10_2);
1911                                const double tmp33 = w8*(-A_11_0 - A_11_2);
1912                                const double tmp34 = w6*(-A_00_0 - A_00_1 - A_00_2 - A_00_3);
1913                                const double tmp35 = w2*(A_01_3 - A_10_1);
1914                                const double tmp36 = w5*(A_01_0 + A_10_3);
1915                                const double tmp37 = w2*(-A_01_3 - A_10_0);
1916                                const double tmp38 = w7*(-A_11_0 - A_11_2);
1917                                const double tmp39 = w5*(-A_01_3 + A_10_1);
1918                                const double tmp40 = w8*(-A_11_1 - A_11_3);
1919                                const double tmp41 = w2*(A_01_0 - A_10_2);
1920                                const double tmp42 = w5*(A_01_1 - A_10_3);
1921                                const double tmp43 = w2*(-A_01_2 + A_10_0);
1922                                const double tmp44 = w5*(A_01_2 - A_10_0);
1923                                const double tmp45 = w2*(-A_01_1 + A_10_3);
1924                                const double tmp46 = w5*(-A_01_0 + A_10_1);
1925                                const double tmp47 = w2*(A_01_3 - A_10_2);
1926                                const double tmp48 = w5*(-A_01_1 - A_10_1);
1927                                const double tmp49 = w2*(A_01_2 + A_10_2);
1928                                const double tmp50 = w5*(-A_01_3 + A_10_2);
1929                                const double tmp51 = w2*(A_01_0 - A_10_1);
1930                                const double tmp52 = w5*(-A_01_2 - A_10_1);
1931                                const double tmp53 = w2*(A_01_1 + A_10_2);
1932                                const double tmp54 = w5*(-A_01_2 - A_10_2);
1933                                const double tmp55 = w2*(A_01_1 + A_10_1);
1934                                EM_S[INDEX2(0,0,4)]+=tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp9;
1935                                EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;
1936                                EM_S[INDEX2(0,2,4)]+=tmp31 + tmp34 + tmp38 + tmp39 + tmp40 + tmp41;
1937                                EM_S[INDEX2(0,3,4)]+=tmp28 + tmp52 + tmp53 + tmp7 + tmp8;
1938                                EM_S[INDEX2(1,0,4)]+=tmp0 + tmp2 + tmp3 + tmp31 + tmp50 + tmp51;
1939                                EM_S[INDEX2(1,1,4)]+=tmp16 + tmp17 + tmp21 + tmp25 + tmp28 + tmp54 + tmp55;
1940                                EM_S[INDEX2(1,2,4)]+=tmp10 + tmp6 + tmp7 + tmp8 + tmp9;
1941                                EM_S[INDEX2(1,3,4)]+=tmp1 + tmp30 + tmp33 + tmp34 + tmp44 + tmp45;
1942                                EM_S[INDEX2(2,0,4)]+=tmp1 + tmp34 + tmp38 + tmp40 + tmp42 + tmp43;
1943                                EM_S[INDEX2(2,1,4)]+=tmp36 + tmp37 + tmp7 + tmp8 + tmp9;
1944                                EM_S[INDEX2(2,2,4)]+=tmp15 + tmp19 + tmp22 + tmp23 + tmp28 + tmp48 + tmp49;
1945                                EM_S[INDEX2(2,3,4)]+=tmp0 + tmp11 + tmp12 + tmp31 + tmp46 + tmp47;
1946                                EM_S[INDEX2(3,0,4)]+=tmp27 + tmp28 + tmp29 + tmp7 + tmp8;
1947                                EM_S[INDEX2(3,1,4)]+=tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35;
1948                                EM_S[INDEX2(3,2,4)]+=tmp0 + tmp1 + tmp11 + tmp12 + tmp13 + tmp14;
1949                                EM_S[INDEX2(3,3,4)]+=tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp9;
1950                            } else { // constant data
1951                                const double A_00 = A_p[INDEX2(0,0,2)];
1952                                const double A_01 = A_p[INDEX2(0,1,2)];
1953                                const double A_10 = A_p[INDEX2(1,0,2)];
1954                                const double A_11 = A_p[INDEX2(1,1,2)];
1955                                const double tmp0 = 6*w1*(A_01 - A_10);
1956                                const double tmp1 = 6*w1*(A_01 + A_10);
1957                                const double tmp2 = 6*w1*(-A_01 - A_10);
1958                                const double tmp3 = 6*w1*(-A_01 + A_10);
1959                                EM_S[INDEX2(0,0,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
1960                                EM_S[INDEX2(0,1,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
1961                                EM_S[INDEX2(0,2,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
1962                                EM_S[INDEX2(0,3,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
1963                                EM_S[INDEX2(1,0,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
1964                                EM_S[INDEX2(1,1,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
1965                                EM_S[INDEX2(1,2,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
1966                                EM_S[INDEX2(1,3,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
1967                                EM_S[INDEX2(2,0,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
1968                                EM_S[INDEX2(2,1,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
1969                                EM_S[INDEX2(2,2,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
1970                                EM_S[INDEX2(2,3,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
1971                                EM_S[INDEX2(3,0,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
1972                                EM_S[INDEX2(3,1,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
1973                                EM_S[INDEX2(3,2,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
1974                                EM_S[INDEX2(3,3,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
1975                            }
1976                        }
1977                        ///////////////
1978                        // process B //
1979                        ///////////////
1980                        if (!B.isEmpty()) {
1981                            add_EM_S=true;
1982                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1983                            if (B.actsExpanded()) {
1984                                const double B_0_0 = B_p[INDEX2(0,0,2)];
1985                                const double B_1_0 = B_p[INDEX2(1,0,2)];
1986                                const double B_0_1 = B_p[INDEX2(0,1,2)];
1987                                const double B_1_1 = B_p[INDEX2(1,1,2)];
1988                                const double B_0_2 = B_p[INDEX2(0,2,2)];
1989                                const double B_1_2 = B_p[INDEX2(1,2,2)];
1990                                const double B_0_3 = B_p[INDEX2(0,3,2)];
1991                                const double B_1_3 = B_p[INDEX2(1,3,2)];
1992                                const double tmp0 = w11*(B_1_0 + B_1_1);
1993                                const double tmp1 = w14*(B_1_2 + B_1_3);
1994                                const double tmp2 = w15*(-B_0_1 - B_0_3);
1995                                const double tmp3 = w10*(-B_0_0 - B_0_2);
1996                                const double tmp4 = w11*(B_1_2 + B_1_3);
1997                                const double tmp5 = w14*(B_1_0 + B_1_1);
1998                                const double tmp6 = w11*(-B_1_2 - B_1_3);
1999                                const double tmp7 = w14*(-B_1_0 - B_1_1);
2000                                const double tmp8 = w11*(-B_1_0 - B_1_1);
2001                                const double tmp9 = w14*(-B_1_2 - B_1_3);
2002                                const double tmp10 = w10*(-B_0_1 - B_0_3);
2003                                const double tmp11 = w15*(-B_0_0 - B_0_2);
2004                                const double tmp12 = w15*(B_0_0 + B_0_2);
2005                                const double tmp13 = w10*(B_0_1 + B_0_3);
2006                                const double tmp14 = w10*(B_0_0 + B_0_2);
2007                                const double tmp15 = w15*(B_0_1 + B_0_3);
2008                                EM_S[INDEX2(0,0,4)]+=B_0_0*w12 + B_0_1*w10 + B_0_2*w15 + B_0_3*w13 + B_1_0*w16 + B_1_1*w14 + B_1_2*w11 + B_1_3*w17;
2009                                EM_S[INDEX2(0,1,4)]+=B_0_0*w10 + B_0_1*w12 + B_0_2*w13 + B_0_3*w15 + tmp0 + tmp1;
2010                                EM_S[INDEX2(0,2,4)]+=B_1_0*w11 + B_1_1*w17 + B_1_2*w16 + B_1_3*w14 + tmp14 + tmp15;
2011                                EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp4 + tmp5;
2012                                EM_S[INDEX2(1,0,4)]+=-B_0_0*w12 - B_0_1*w10 - B_0_2*w15 - B_0_3*w13 + tmp0 + tmp1;
2013                                EM_S[INDEX2(1,1,4)]+=-B_0_0*w10 - B_0_1*w12 - B_0_2*w13 - B_0_3*w15 + B_1_0*w14 + B_1_1*w16 + B_1_2*w17 + B_1_3*w11;
2014                                EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2015                                EM_S[INDEX2(1,3,4)]+=B_1_0*w17 + B_1_1*w11 + B_1_2*w14 + B_1_3*w16 + tmp10 + tmp11;
2016                                EM_S[INDEX2(2,0,4)]+=-B_1_0*w16 - B_1_1*w14 - B_1_2*w11 - B_1_3*w17 + tmp14 + tmp15;
2017                                EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2018                                EM_S[INDEX2(2,2,4)]+=B_0_0*w15 + B_0_1*w13 + B_0_2*w12 + B_0_3*w10 - B_1_0*w11 - B_1_1*w17 - B_1_2*w16 - B_1_3*w14;
2019                                EM_S[INDEX2(2,3,4)]+=B_0_0*w13 + B_0_1*w15 + B_0_2*w10 + B_0_3*w12 + tmp6 + tmp7;
2020                                EM_S[INDEX2(3,0,4)]+=tmp2 + tmp3 + tmp8 + tmp9;
2021                                EM_S[INDEX2(3,1,4)]+=-B_1_0*w14 - B_1_1*w16 - B_1_2*w17 - B_1_3*w11 + tmp10 + tmp11;
2022                                EM_S[INDEX2(3,2,4)]+=-B_0_0*w15 - B_0_1*w13 - B_0_2*w12 - B_0_3*w10 + tmp6 + tmp7;
2023                                EM_S[INDEX2(3,3,4)]+=-B_0_0*w13 - B_0_1*w15 - B_0_2*w10 - B_0_3*w12 - B_1_0*w17 - B_1_1*w11 - B_1_2*w14 - B_1_3*w16;
2024                            } else { // constant data
2025                                const double B_0 = B_p[0];
2026                                const double B_1 = B_p[1];
2027                                EM_S[INDEX2(0,0,4)]+= 2*B_0*w18 + 2*B_1*w19;
2028                                EM_S[INDEX2(0,1,4)]+= 2*B_0*w18 +   B_1*w19;
2029                                EM_S[INDEX2(0,2,4)]+=   B_0*w18 + 2*B_1*w19;
2030                                EM_S[INDEX2(0,3,4)]+=   B_0*w18 +   B_1*w19;
2031                                EM_S[INDEX2(1,0,4)]+=-2*B_0*w18 +   B_1*w19;
2032                                EM_S[INDEX2(1,1,4)]+=-2*B_0*w18 + 2*B_1*w19;
2033                                EM_S[INDEX2(1,2,4)]+=  -B_0*w18 +   B_1*w19;
2034                                EM_S[INDEX2(1,3,4)]+=  -B_0*w18 + 2*B_1*w19;
2035                                EM_S[INDEX2(2,0,4)]+=   B_0*w18 - 2*B_1*w19;
2036                                EM_S[INDEX2(2,1,4)]+=   B_0*w18 -   B_1*w19;
2037                                EM_S[INDEX2(2,2,4)]+= 2*B_0*w18 - 2*B_1*w19;
2038                                EM_S[INDEX2(2,3,4)]+= 2*B_0*w18 -   B_1*w19;
2039                                EM_S[INDEX2(3,0,4)]+=  -B_0*w18 -   B_1*w19;
2040                                EM_S[INDEX2(3,1,4)]+=  -B_0*w18 - 2*B_1*w19;
2041                                EM_S[INDEX2(3,2,4)]+=-2*B_0*w18 -   B_1*w19;
2042                                EM_S[INDEX2(3,3,4)]+=-2*B_0*w18 - 2*B_1*w19;
2043                            }
2044                        }
2045                        ///////////////
2046                        // process C //
2047                        ///////////////
2048                        if (!C.isEmpty()) {
2049                            add_EM_S=true;
2050                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2051                            if (C.actsExpanded()) {
2052                                const double C_0_0 = C_p[INDEX2(0,0,2)];
2053                                const double C_1_0 = C_p[INDEX2(1,0,2)];
2054                                const double C_0_1 = C_p[INDEX2(0,1,2)];
2055                                const double C_1_1 = C_p[INDEX2(1,1,2)];
2056                                const double C_0_2 = C_p[INDEX2(0,2,2)];
2057                                const double C_1_2 = C_p[INDEX2(1,2,2)];
2058                                const double C_0_3 = C_p[INDEX2(0,3,2)];
2059                                const double C_1_3 = C_p[INDEX2(1,3,2)];
2060                                const double tmp0 = w11*(C_1_0 + C_1_1);
2061                                const double tmp1 = w14*(C_1_2 + C_1_3);
2062                                const double tmp2 = w15*(C_0_0 + C_0_2);
2063                                const double tmp3 = w10*(C_0_1 + C_0_3);
2064                                const double tmp4 = w11*(-C_1_0 - C_1_1);
2065                                const double tmp5 = w14*(-C_1_2 - C_1_3);
2066                                const double tmp6 = w11*(-C_1_2 - C_1_3);
2067                                const double tmp7 = w14*(-C_1_0 - C_1_1);
2068                                const double tmp8 = w11*(C_1_2 + C_1_3);
2069                                const double tmp9 = w14*(C_1_0 + C_1_1);
2070                                const double tmp10 = w10*(-C_0_1 - C_0_3);
2071                                const double tmp11 = w15*(-C_0_0 - C_0_2);
2072                                const double tmp12 = w15*(-C_0_1 - C_0_3);
2073                                const double tmp13 = w10*(-C_0_0 - C_0_2);
2074                                const double tmp14 = w10*(C_0_0 + C_0_2);
2075                                const double tmp15 = w15*(C_0_1 + C_0_3);
2076                                EM_S[INDEX2(0,0,4)]+=C_0_0*w12 + C_0_1*w10 + C_0_2*w15 + C_0_3*w13 + C_1_0*w16 + C_1_1*w14 + C_1_2*w11 + C_1_3*w17;
2077                                EM_S[INDEX2(0,1,4)]+=-C_0_0*w12 - C_0_1*w10 - C_0_2*w15 - C_0_3*w13 + tmp0 + tmp1;
2078                                EM_S[INDEX2(0,2,4)]+=-C_1_0*w16 - C_1_1*w14 - C_1_2*w11 - C_1_3*w17 + tmp14 + tmp15;
2079                                EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp4 + tmp5;
2080                                EM_S[INDEX2(1,0,4)]+=C_0_0*w10 + C_0_1*w12 + C_0_2*w13 + C_0_3*w15 + tmp0 + tmp1;
2081                                EM_S[INDEX2(1,1,4)]+=-C_0_0*w10 - C_0_1*w12 - C_0_2*w13 - C_0_3*w15 + C_1_0*w14 + C_1_1*w16 + C_1_2*w17 + C_1_3*w11;
2082                                EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2083                                EM_S[INDEX2(1,3,4)]+=-C_1_0*w14 - C_1_1*w16 - C_1_2*w17 - C_1_3*w11 + tmp10 + tmp11;
2084                                EM_S[INDEX2(2,0,4)]+=C_1_0*w11 + C_1_1*w17 + C_1_2*w16 + C_1_3*w14 + tmp14 + tmp15;
2085                                EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2086                                EM_S[INDEX2(2,2,4)]+=C_0_0*w15 + C_0_1*w13 + C_0_2*w12 + C_0_3*w10 - C_1_0*w11 - C_1_1*w17 - C_1_2*w16 - C_1_3*w14;
2087                                EM_S[INDEX2(2,3,4)]+=-C_0_0*w15 - C_0_1*w13 - C_0_2*w12 - C_0_3*w10 + tmp6 + tmp7;
2088                                EM_S[INDEX2(3,0,4)]+=tmp2 + tmp3 + tmp8 + tmp9;
2089                                EM_S[INDEX2(3,1,4)]+=C_1_0*w17 + C_1_1*w11 + C_1_2*w14 + C_1_3*w16 + tmp10 + tmp11;
2090                                EM_S[INDEX2(3,2,4)]+=C_0_0*w13 + C_0_1*w15 + C_0_2*w10 + C_0_3*w12 + tmp6 + tmp7;
2091                                EM_S[INDEX2(3,3,4)]+=-C_0_0*w13 - C_0_1*w15 - C_0_2*w10 - C_0_3*w12 - C_1_0*w17 - C_1_1*w11 - C_1_2*w14 - C_1_3*w16;
2092                            } else { // constant data
2093                                const double C_0 = C_p[0];
2094                                const double C_1 = C_p[1];
2095                                EM_S[INDEX2(0,0,4)]+= 2*C_0*w18 + 2*C_1*w19;
2096                                EM_S[INDEX2(0,1,4)]+=-2*C_0*w18 +   C_1*w19;
2097                                EM_S[INDEX2(0,2,4)]+=   C_0*w18 - 2*C_1*w19;
2098                                EM_S[INDEX2(0,3,4)]+=  -C_0*w18 -   C_1*w19;
2099                                EM_S[INDEX2(1,0,4)]+= 2*C_0*w18 +   C_1*w19;
2100                                EM_S[INDEX2(1,1,4)]+=-2*C_0*w18 + 2*C_1*w19;
2101                                EM_S[INDEX2(1,2,4)]+=   C_0*w18 -   C_1*w19;
2102                                EM_S[INDEX2(1,3,4)]+=  -C_0*w18 - 2*C_1*w19;
2103                                EM_S[INDEX2(2,0,4)]+=   C_0*w18 + 2*C_1*w19;
2104                                EM_S[INDEX2(2,1,4)]+=  -C_0*w18 +   C_1*w19;
2105                                EM_S[INDEX2(2,2,4)]+= 2*C_0*w18 - 2*C_1*w19;
2106                                EM_S[INDEX2(2,3,4)]+=-2*C_0*w18 -   C_1*w19;
2107                                EM_S[INDEX2(3,0,4)]+=   C_0*w18 +   C_1*w19;
2108                                EM_S[INDEX2(3,1,4)]+=  -C_0*w18 + 2*C_1*w19;
2109                                EM_S[INDEX2(3,2,4)]+= 2*C_0*w18 -   C_1*w19;
2110                                EM_S[INDEX2(3,3,4)]+=-2*C_0*w18 - 2*C_1*w19;
2111                            }
2112                        }
2113                        ///////////////
2114                        // process D //
2115                        ///////////////
2116                        if (!D.isEmpty()) {
2117                            add_EM_S=true;
2118                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2119                            if (D.actsExpanded()) {
2120                                const double D_0 = D_p[0];
2121                                const double D_1 = D_p[1];
2122                                const double D_2 = D_p[2];
2123                                const double D_3 = D_p[3];
2124                                const double tmp0 = w21*(D_2 + D_3);
2125                                const double tmp1 = w20*(D_0 + D_1);
2126                                const double tmp2 = w22*(D_0 + D_1 + D_2 + D_3);
2127                                const double tmp3 = w21*(D_0 + D_1);
2128                                const double tmp4 = w20*(D_2 + D_3);
2129                                const double tmp5 = w22*(D_1 + D_2);
2130                                const double tmp6 = w21*(D_0 + D_2);
2131                                const double tmp7 = w20*(D_1 + D_3);
2132                                const double tmp8 = w21*(D_1 + D_3);
2133                                const double tmp9 = w20*(D_0 + D_2);
2134                                const double tmp10 = w22*(D_0 + D_3);
2135                                EM_S[INDEX2(0,0,4)]+=D_0*w23 + D_3*w24 + tmp5;
2136                                EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1;
2137                                EM_S[INDEX2(0,2,4)]+=tmp8 + tmp9;
2138                                EM_S[INDEX2(0,3,4)]+=tmp2;
2139                                EM_S[INDEX2(1,0,4)]+=tmp0 + tmp1;
2140                                EM_S[INDEX2(1,1,4)]+=D_1*w23 + D_2*w24 + tmp10;
2141                                EM_S[INDEX2(1,2,4)]+=tmp2;
2142                                EM_S[INDEX2(1,3,4)]+=tmp6 + tmp7;
2143                                EM_S[INDEX2(2,0,4)]+=tmp8 + tmp9;
2144                                EM_S[INDEX2(2,1,4)]+=tmp2;
2145                                EM_S[INDEX2(2,2,4)]+=D_1*w24 + D_2*w23 + tmp10;
2146                                EM_S[INDEX2(2,3,4)]+=tmp3 + tmp4;
2147                                EM_S[INDEX2(3,0,4)]+=tmp2;
2148                                EM_S[INDEX2(3,1,4)]+=tmp6 + tmp7;
2149                                EM_S[INDEX2(3,2,4)]+=tmp3 + tmp4;
2150                                EM_S[INDEX2(3,3,4)]+=D_0*w24 + D_3*w23 + tmp5;
2151                            } else { // constant data
2152                                const double D_0 = D_p[0];
2153                                EM_S[INDEX2(0,0,4)]+=16*D_0*w22;
2154                                EM_S[INDEX2(0,1,4)]+=8*D_0*w22;
2155                                EM_S[INDEX2(0,2,4)]+=8*D_0*w22;
2156                                EM_S[INDEX2(0,3,4)]+=4*D_0*w22;
2157                                EM_S[INDEX2(1,0,4)]+=8*D_0*w22;
2158                                EM_S[INDEX2(1,1,4)]+=16*D_0*w22;
2159                                EM_S[INDEX2(1,2,4)]+=4*D_0*w22;
2160                                EM_S[INDEX2(1,3,4)]+=8*D_0*w22;
2161                                EM_S[INDEX2(2,0,4)]+=8*D_0*w22;
2162                                EM_S[INDEX2(2,1,4)]+=4*D_0*w22;
2163                                EM_S[INDEX2(2,2,4)]+=16*D_0*w22;
2164                                EM_S[INDEX2(2,3,4)]+=8*D_0*w22;
2165                                EM_S[INDEX2(3,0,4)]+=4*D_0*w22;
2166                                EM_S[INDEX2(3,1,4)]+=8*D_0*w22;
2167                                EM_S[INDEX2(3,2,4)]+=8*D_0*w22;
2168                                EM_S[INDEX2(3,3,4)]+=16*D_0*w22;
2169                            }
2170                        }
2171                        ///////////////
2172                        // process X //
2173                        ///////////////
2174                        if (!X.isEmpty()) {
2175                            add_EM_F=true;
2176                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2177                            if (X.actsExpanded()) {
2178                                const double X_0_0 = X_p[INDEX2(0,0,2)];
2179                                const double X_1_0 = X_p[INDEX2(1,0,2)];
2180                                const double X_0_1 = X_p[INDEX2(0,1,2)];
2181                                const double X_1_1 = X_p[INDEX2(1,1,2)];
2182                                const double X_0_2 = X_p[INDEX2(0,2,2)];
2183                                const double X_1_2 = X_p[INDEX2(1,2,2)];
2184                                const double X_0_3 = X_p[INDEX2(0,3,2)];
2185                                const double X_1_3 = X_p[INDEX2(1,3,2)];
2186                                const double tmp0 = 6*w15*(X_0_2 + X_0_3);
2187                                const double tmp1 = 6*w10*(X_0_0 + X_0_1);
2188                                const double tmp2 = 6*w11*(X_1_0 + X_1_2);
2189                                const double tmp3 = 6*w14*(X_1_1 + X_1_3);
2190                                const double tmp4 = 6*w11*(X_1_1 + X_1_3);
2191                                const double tmp5 = w25*(X_0_0 + X_0_1);
2192                                const double tmp6 = w26*(X_0_2 + X_0_3);
2193                                const double tmp7 = 6*w14*(X_1_0 + X_1_2);
2194                                const double tmp8 = w27*(X_1_0 + X_1_2);
2195                                const double tmp9 = w28*(X_1_1 + X_1_3);
2196                                const double tmp10 = w25*(-X_0_2 - X_0_3);
2197                                const double tmp11 = w26*(-X_0_0 - X_0_1);
2198                                const double tmp12 = w27*(X_1_1 + X_1_3);
2199                                const double tmp13 = w28*(X_1_0 + X_1_2);
2200                                const double tmp14 = w25*(X_0_2 + X_0_3);
2201                                const double tmp15 = w26*(X_0_0 + X_0_1);
2202                                EM_F[0]+=tmp0 + tmp1 + tmp2 + tmp3;
2203                                EM_F[1]+=tmp4 + tmp5 + tmp6 + tmp7;
2204                                EM_F[2]+=tmp10 + tmp11 + tmp8 + tmp9;
2205                                EM_F[3]+=tmp12 + tmp13 + tmp14 + tmp15;
2206                            } else { // constant data
2207                                const double X_0 = X_p[0];
2208                                const double X_1 = X_p[1];
2209                                EM_F[0]+= 6*X_0*w18 + 6*X_1*w19;
2210                                EM_F[1]+=-6*X_0*w18 + 6*X_1*w19;
2211                                EM_F[2]+= 6*X_0*w18 - 6*X_1*w19;
2212                                EM_F[3]+=-6*X_0*w18 - 6*X_1*w19;
2213                            }
2214                        }
2215                        ///////////////
2216                        // process Y //
2217                        ///////////////
2218                        if (!Y.isEmpty()) {
2219                            add_EM_F=true;
2220                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2221                            if (Y.actsExpanded()) {
2222                                const double Y_0 = Y_p[0];
2223                                const double Y_1 = Y_p[1];
2224                                const double Y_2 = Y_p[2];
2225                                const double Y_3 = Y_p[3];
2226                                const double tmp0 = 6*w22*(Y_1 + Y_2);
2227                                const double tmp1 = 6*w22*(Y_0 + Y_3);
2228                                EM_F[0]+=6*Y_0*w20 + 6*Y_3*w21 + tmp0;
2229                                EM_F[1]+=6*Y_1*w20 + 6*Y_2*w21 + tmp1;
2230                                EM_F[2]+=6*Y_1*w21 + 6*Y_2*w20 + tmp1;
2231                                EM_F[3]+=6*Y_0*w21 + 6*Y_3*w20 + tmp0;
2232                            } else { // constant data
2233                                EM_F[0]+=36*Y_p[0]*w22;
2234                                EM_F[1]+=36*Y_p[0]*w22;
2235                                EM_F[2]+=36*Y_p[0]*w22;
2236                                EM_F[3]+=36*Y_p[0]*w22;
2237                            }
2238                        }
2239    
2240                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2241                        const index_t firstNode=m_NN[0]*k1+k0;
2242                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2243                    } // end k0 loop
2244                } // end k1 loop
2245            } // end of colouring
2246        } // end of parallel region
2247    }
2248    
2249    //protected
2250    void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2251            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2252            const escript::Data& C, const escript::Data& D,
2253            const escript::Data& X, const escript::Data& Y) const
2254    {
2255        const double w0 = 1./4;
2256        const double w1 = m_dx[0]/8;
2257        const double w2 = m_dx[1]/8;
2258        const double w3 = m_dx[0]*m_dx[1]/16;
2259        const double w4 = m_dx[0]/(4*m_dx[1]);
2260        const double w5 = m_dx[1]/(4*m_dx[0]);
2261    
2262        rhs.requireWrite();
2263    #pragma omp parallel
2264        {
2265            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2266    #pragma omp for
2267                for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
2268                    for (index_t k0=0; k0<m_NE[0]; ++k0)  {
2269                        bool add_EM_S=false;
2270                        bool add_EM_F=false;
2271                        vector<double> EM_S(4*4, 0);
2272                        vector<double> EM_F(4, 0);
2273                        const index_t e = k0 + m_NE[0]*k1;
2274                        ///////////////
2275                        // process A //
2276                        ///////////////
2277                        if (!A.isEmpty()) {
2278                            add_EM_S=true;
2279                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2280                            const double A_00 = A_p[INDEX2(0,0,2)];
2281                            const double A_10 = A_p[INDEX2(1,0,2)];
2282                            const double A_01 = A_p[INDEX2(0,1,2)];
2283                            const double A_11 = A_p[INDEX2(1,1,2)];
2284                            const double tmp0 = (A_01 + A_10)*w0;
2285                            const double tmp1 = A_00*w5;
2286                            const double tmp2 = A_01*w0;
2287                            const double tmp3 = A_10*w0;
2288                            const double tmp4 = A_11*w4;
2289                            EM_S[INDEX2(0,0,4)]+=tmp4 + tmp0 + tmp1;
2290                            EM_S[INDEX2(1,0,4)]+=tmp4 - tmp1 + tmp3 - tmp2;
2291                            EM_S[INDEX2(2,0,4)]+=tmp2 - tmp3 - tmp4 + tmp1;
2292                            EM_S[INDEX2(3,0,4)]+=-tmp1 - tmp4 - tmp0;
2293                            EM_S[INDEX2(0,1,4)]+=tmp4 - tmp1 + tmp2 - tmp3;
2294                            EM_S[INDEX2(1,1,4)]+=tmp4 + tmp1 - tmp0;
2295                            EM_S[INDEX2(2,1,4)]+=-tmp1 + tmp0 - tmp4;
2296                            EM_S[INDEX2(3,1,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;
2297                            EM_S[INDEX2(0,2,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;
2298                            EM_S[INDEX2(1,2,4)]+=-tmp1 + tmp0 - tmp4;
2299                            EM_S[INDEX2(2,2,4)]+=tmp4 + tmp1 - tmp0;
2300                            EM_S[INDEX2(3,2,4)]+=tmp4 - tmp1 + tmp2 - tmp3;
2301                            EM_S[INDEX2(0,3,4)]+=-tmp1 - tmp4 - tmp0;
2302                            EM_S[INDEX2(1,3,4)]+=tmp2 - tmp3 - tmp4 + tmp1;
2303                            EM_S[INDEX2(2,3,4)]+=tmp4 - tmp1 + tmp3 - tmp2;
2304                            EM_S[INDEX2(3,3,4)]+=tmp4 + tmp0 + tmp1;
2305                        }
2306                        ///////////////
2307                        // process B //
2308                        ///////////////
2309                        if (!B.isEmpty()) {
2310                            add_EM_S=true;
2311                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2312                            const double tmp0 = B_p[0]*w2;
2313                            const double tmp1 = B_p[1]*w1;
2314                            EM_S[INDEX2(0,0,4)]+=-tmp0 - tmp1;
2315                            EM_S[INDEX2(1,0,4)]+= tmp0 - tmp1;
2316                            EM_S[INDEX2(2,0,4)]+= tmp1 - tmp0;
2317                            EM_S[INDEX2(3,0,4)]+= tmp0 + tmp1;
2318                            EM_S[INDEX2(0,1,4)]+=-tmp0 - tmp1;
2319                            EM_S[INDEX2(1,1,4)]+= tmp0 - tmp1;
2320                            EM_S[INDEX2(2,1,4)]+= tmp1 - tmp0;
2321                            EM_S[INDEX2(3,1,4)]+= tmp0 + tmp1;
2322                            EM_S[INDEX2(0,2,4)]+=-tmp0 - tmp1;
2323                            EM_S[INDEX2(1,2,4)]+= tmp0 - tmp1;
2324                            EM_S[INDEX2(2,2,4)]+= tmp1 - tmp0;
2325                            EM_S[INDEX2(3,2,4)]+= tmp0 + tmp1;
2326                            EM_S[INDEX2(0,3,4)]+=-tmp0 - tmp1;
2327                            EM_S[INDEX2(1,3,4)]+= tmp0 - tmp1;
2328                            EM_S[INDEX2(2,3,4)]+= tmp1 - tmp0;
2329                            EM_S[INDEX2(3,3,4)]+= tmp0 + tmp1;
2330                        }
2331                        ///////////////
2332                        // process C //
2333                        ///////////////
2334                        if (!C.isEmpty()) {
2335                            add_EM_S=true;
2336                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2337                            const double tmp0 = C_p[0]*w2;
2338                            const double tmp1 = C_p[1]*w1;
2339                            EM_S[INDEX2(0,0,4)]+=-tmp1 - tmp0;
2340                            EM_S[INDEX2(1,0,4)]+=-tmp1 - tmp0;
2341                            EM_S[INDEX2(2,0,4)]+=-tmp1 - tmp0;
2342                            EM_S[INDEX2(3,0,4)]+=-tmp1 - tmp0;
2343                            EM_S[INDEX2(0,1,4)]+= tmp0 - tmp1;
2344                            EM_S[INDEX2(1,1,4)]+= tmp0 - tmp1;
2345                            EM_S[INDEX2(2,1,4)]+= tmp0 - tmp1;
2346                            EM_S[INDEX2(3,1,4)]+= tmp0 - tmp1;
2347                            EM_S[INDEX2(0,2,4)]+= tmp1 - tmp0;
2348                            EM_S[INDEX2(1,2,4)]+= tmp1 - tmp0;
2349                            EM_S[INDEX2(2,2,4)]+= tmp1 - tmp0;
2350                            EM_S[INDEX2(3,2,4)]+= tmp1 - tmp0;
2351                            EM_S[INDEX2(0,3,4)]+= tmp0 + tmp1;
2352                            EM_S[INDEX2(1,3,4)]+= tmp0 + tmp1;
2353                            EM_S[INDEX2(2,3,4)]+= tmp0 + tmp1;
2354                            EM_S[INDEX2(3,3,4)]+= tmp0 + tmp1;
2355                        }
2356                        ///////////////
2357                        // process D //
2358                        ///////////////
2359                        if (!D.isEmpty()) {
2360                            add_EM_S=true;
2361                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2362                            EM_S[INDEX2(0,0,4)]+=D_p[0]*w3;
2363                            EM_S[INDEX2(1,0,4)]+=D_p[0]*w3;
2364                            EM_S[INDEX2(2,0,4)]+=D_p[0]*w3;
2365                            EM_S[INDEX2(3,0,4)]+=D_p[0]*w3;
2366                            EM_S[INDEX2(0,1,4)]+=D_p[0]*w3;
2367                            EM_S[INDEX2(1,1,4)]+=D_p[0]*w3;
2368                            EM_S[INDEX2(2,1,4)]+=D_p[0]*w3;
2369                            EM_S[INDEX2(3,1,4)]+=D_p[0]*w3;
2370                            EM_S[INDEX2(0,2,4)]+=D_p[0]*w3;
2371                            EM_S[INDEX2(1,2,4)]+=D_p[0]*w3;
2372                            EM_S[INDEX2(2,2,4)]+=D_p[0]*w3;
2373                            EM_S[INDEX2(3,2,4)]+=D_p[0]*w3;
2374                            EM_S[INDEX2(0,3,4)]+=D_p[0]*w3;
2375                            EM_S[INDEX2(1,3,4)]+=D_p[0]*w3;
2376                            EM_S[INDEX2(2,3,4)]+=D_p[0]*w3;
2377                            EM_S[INDEX2(3,3,4)]+=D_p[0]*w3;
2378                        }
2379                        ///////////////
2380                        // process X //
2381                        ///////////////
2382                        if (!X.isEmpty()) {
2383                            add_EM_F=true;
2384                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2385                            const double wX0 = 4*X_p[0]*w2;
2386                            const double wX1 = 4*X_p[1]*w1;
2387                            EM_F[0]+=-wX0 - wX1;
2388                            EM_F[1]+=-wX1 + wX0;
2389                            EM_F[2]+=-wX0 + wX1;
2390                            EM_F[3]+= wX0 + wX1;
2391                        }
2392                        ///////////////
2393                        // process Y //
2394                        ///////////////
2395                        if (!Y.isEmpty()) {
2396                            add_EM_F=true;
2397                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2398                            EM_F[0]+=4*Y_p[0]*w3;
2399                            EM_F[1]+=4*Y_p[0]*w3;
2400                            EM_F[2]+=4*Y_p[0]*w3;
2401                            EM_F[3]+=4*Y_p[0]*w3;
2402                        }
2403    
2404                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2405                        const index_t firstNode=m_NN[0]*k1+k0;
2406                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2407                    } // end k0 loop
2408                } // end k1 loop
2409            } // end of colouring
2410        } // end of parallel region
2411    }
2412    
2413    //protected
2414    void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2415            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2416            const escript::Data& C, const escript::Data& D,
2417            const escript::Data& X, const escript::Data& Y) const
2418    {
2419        dim_t numEq, numComp;
2420        if (!mat)
2421            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
2422        else {
2423            numEq=mat->logical_row_block_size;
2424            numComp=mat->logical_col_block_size;
2425        }
2426        const double SQRT3 = 1.73205080756887719318;
2427        const double w1 = 1.0/24;
2428        const double w5 = -SQRT3/24 + 1.0/12;
2429        const double w2 = -SQRT3/24 - 1.0/12;
2430        const double w19 = -m_dx[0]/12;
2431        const double w11 = w19*(SQRT3 + 3)/12;
2432        const double w14 = w19*(-SQRT3 + 3)/12;
2433        const double w16 = w19*(5*SQRT3 + 9)/12;
2434        const double w17 = w19*(-5*SQRT3 + 9)/12;
2435        const double w27 = w19*(-SQRT3 - 3)/2;
2436        const double w28 = w19*(SQRT3 - 3)/2;
2437        const double w18 = -m_dx[1]/12;
2438        const double w10 = w18*(SQRT3 + 3)/12;
2439        const double w15 = w18*(-SQRT3 + 3)/12;
2440        const double w12 = w18*(5*SQRT3 + 9)/12;
2441        const double w13 = w18*(-5*SQRT3 + 9)/12;
2442        const double w25 = w18*(-SQRT3 - 3)/2;
2443        const double w26 = w18*(SQRT3 - 3)/2;
2444        const double w22 = m_dx[0]*m_dx[1]/144;
2445        const double w20 = w22*(SQRT3 + 2);
2446        const double w21 = w22*(-SQRT3 + 2);
2447        const double w23 = w22*(4*SQRT3 + 7);
2448        const double w24 = w22*(-4*SQRT3 + 7);
2449        const double w3 = m_dx[0]/(24*m_dx[1]);
2450        const double w7 = w3*(SQRT3 + 2);
2451        const double w8 = w3*(-SQRT3 + 2);
2452        const double w6 = -m_dx[1]/(24*m_dx[0]);
2453        const double w0 = w6*(SQRT3 + 2);
2454        const double w4 = w6*(-SQRT3 + 2);
2455    
2456        rhs.requireWrite();
2457    #pragma omp parallel
2458        {
2459            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2460    #pragma omp for
2461                for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
2462                    for (index_t k0=0; k0<m_NE[0]; ++k0)  {
2463                        bool add_EM_S=false;
2464                        bool add_EM_F=false;
2465                        vector<double> EM_S(4*4*numEq*numComp, 0);
2466                        vector<double> EM_F(4*numEq, 0);
2467                        const index_t e = k0 + m_NE[0]*k1;
2468                        ///////////////
2469                        // process A //
2470                        ///////////////
2471                        if (!A.isEmpty()) {
2472                            add_EM_S = true;
2473                            const double* A_p = const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2474                            if (A.actsExpanded()) {
2475                                for (index_t k=0; k<numEq; k++) {
2476                                    for (index_t m=0; m<numComp; m++) {
2477                                        const double A_00_0 = A_p[INDEX5(k,0,m,0,0,numEq,2,numComp,2)];
2478                                        const double A_01_0 = A_p[INDEX5(k,0,m,1,0,numEq,2,numComp,2)];
2479                                        const double A_10_0 = A_p[INDEX5(k,1,m,0,0,numEq,2,numComp,2)];
2480                                        const double A_11_0 = A_p[INDEX5(k,1,m,1,0,numEq,2,numComp,2)];
2481                                        const double A_00_1 = A_p[INDEX5(k,0,m,0,1,numEq,2,numComp,2)];
2482                                        const double A_01_1 = A_p[INDEX5(k,0,m,1,1,numEq,2,numComp,2)];
2483                                        const double A_10_1 = A_p[INDEX5(k,1,m,0,1,numEq,2,numComp,2)];
2484                                        const double A_11_1 = A_p[INDEX5(k,1,m,1,1,numEq,2,numComp,2)];
2485                                        const double A_00_2 = A_p[INDEX5(k,0,m,0,2,numEq,2,numComp,2)];
2486                                        const double A_01_2 = A_p[INDEX5(k,0,m,1,2,numEq,2,numComp,2)];
2487                                        const double A_10_2 = A_p[INDEX5(k,1,m,0,2,numEq,2,numComp,2)];
2488                                        const double A_11_2 = A_p[INDEX5(k,1,m,1,2,numEq,2,numComp,2)];
2489                                        const double A_00_3 = A_p[INDEX5(k,0,m,0,3,numEq,2,numComp,2)];
2490                                        const double A_01_3 = A_p[INDEX5(k,0,m,1,3,numEq,2,numComp,2)];
2491                                        const double A_10_3 = A_p[INDEX5(k,1,m,0,3,numEq,2,numComp,2)];
2492                                        const double A_11_3 = A_p[INDEX5(k,1,m,1,3,numEq,2,numComp,2)];
2493                                        const double tmp0 = w3*(A_11_0 + A_11_1 + A_11_2 + A_11_3);
2494                                        const double tmp1 = w1*(A_01_0 + A_01_3 - A_10_1 - A_10_2);
2495                                        const double tmp2 = w4*(A_00_2 + A_00_3);
2496                                        const double tmp3 = w0*(A_00_0 + A_00_1);
2497                                        const double tmp4 = w5*(A_01_2 - A_10_3);
2498                                        const double tmp5 = w2*(-A_01_1 + A_10_0);
2499                                        const double tmp6 = w5*(A_01_3 + A_10_0);
2500                                        const double tmp7 = w3*(-A_11_0 - A_11_1 - A_11_2 - A_11_3);
2501                                        const double tmp8 = w6*(A_00_0 + A_00_1 + A_00_2 + A_00_3);
2502                                        const double tmp9 = w1*(A_01_1 + A_01_2 + A_10_1 + A_10_2);
2503                                        const double tmp10 = w2*(-A_01_0 - A_10_3);
2504                                        const double tmp11 = w4*(A_00_0 + A_00_1);
2505                                        const double tmp12 = w0*(A_00_2 + A_00_3);
2506                                        const double tmp13 = w5*(A_01_1 - A_10_0);
2507                                        const double tmp14 = w2*(-A_01_2 + A_10_3);
2508                                        const double tmp15 = w7*(A_11_0 + A_11_2);
2509                                        const double tmp16 = w4*(-A_00_2 - A_00_3);
2510                                        const double tmp17 = w0*(-A_00_0 - A_00_1);
2511                                        const double tmp18 = w5*(A_01_3 + A_10_3);
2512                                        const double tmp19 = w8*(A_11_1 + A_11_3);
2513                                        const double tmp20 = w2*(-A_01_0 - A_10_0);
2514                                        const double tmp21 = w7*(A_11_1 + A_11_3);
2515                                        const double tmp22 = w4*(-A_00_0 - A_00_1);
2516                                        const double tmp23 = w0*(-A_00_2 - A_00_3);
2517                                        const double tmp24 = w5*(A_01_0 + A_10_0);
2518                                        const double tmp25 = w8*(A_11_0 + A_11_2);
2519                                        const double tmp26 = w2*(-A_01_3 - A_10_3);
2520                                        const double tmp27 = w5*(-A_01_1 - A_10_2);
2521                                        const double tmp28 = w1*(-A_01_0 - A_01_3 - A_10_0 - A_10_3);
2522                                        const double tmp29 = w2*(A_01_2 + A_10_1);
2523                                        const double tmp30 = w7*(-A_11_1 - A_11_3);
2524                                        const double tmp31 = w1*(-A_01_1 - A_01_2 + A_10_0 + A_10_3);
2525                                        const double tmp32 = w5*(-A_01_0 + A_10_2);
2526                                        const double tmp33 = w8*(-A_11_0 - A_11_2);
2527                                        const double tmp34 = w6*(-A_00_0 - A_00_1 - A_00_2 - A_00_3);
2528                                        const double tmp35 = w2*(A_01_3 - A_10_1);
2529                                        const double tmp36 = w5*(A_01_0 + A_10_3);
2530                                        const double tmp37 = w2*(-A_01_3 - A_10_0);
2531                                        const double tmp38 = w7*(-A_11_0 - A_11_2);
2532                                        const double tmp39 = w5*(-A_01_3 + A_10_1);
2533                                        const double tmp40 = w8*(-A_11_1 - A_11_3);
2534                                        const double tmp41 = w2*(A_01_0 - A_10_2);
2535                                        const double tmp42 = w5*(A_01_1 - A_10_3);
2536                                        const double tmp43 = w2*(-A_01_2 + A_10_0);
2537                                        const double tmp44 = w5*(A_01_2 - A_10_0);
2538                                        const double tmp45 = w2*(-A_01_1 + A_10_3);
2539                                        const double tmp46 = w5*(-A_01_0 + A_10_1);
2540                                        const double tmp47 = w2*(A_01_3 - A_10_2);
2541                                        const double tmp48 = w5*(-A_01_1 - A_10_1);
2542                                        const double tmp49 = w2*(A_01_2 + A_10_2);
2543                                        const double tmp50 = w5*(-A_01_3 + A_10_2);
2544                                        const double tmp51 = w2*(A_01_0 - A_10_1);
2545                                        const double tmp52 = w5*(-A_01_2 - A_10_1);
2546                                        const double tmp53 = w2*(A_01_1 + A_10_2);
2547                                        const double tmp54 = w5*(-A_01_2 - A_10_2);
2548                                        const double tmp55 = w2*(A_01_1 + A_10_1);
2549                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp9;
2550                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;
2551                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp31 + tmp34 + tmp38 + tmp39 + tmp40 + tmp41;
2552                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp28 + tmp52 + tmp53 + tmp7 + tmp8;
2553                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp2 + tmp3 + tmp31 + tmp50 + tmp51;
2554                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp16 + tmp17 + tmp21 + tmp25 + tmp28 + tmp54 + tmp55;
2555                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10 + tmp6 + tmp7 + tmp8 + tmp9;
2556                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1 + tmp30 + tmp33 + tmp34 + tmp44 + tmp45;
2557                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1 + tmp34 + tmp38 + tmp40 + tmp42 + tmp43;
2558                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp36 + tmp37 + tmp7 + tmp8 + tmp9;
2559                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp15 + tmp19 + tmp22 + tmp23 + tmp28 + tmp48 + tmp49;
2560                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0 + tmp11 + tmp12 + tmp31 + tmp46 + tmp47;
2561                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp27 + tmp28 + tmp29 + tmp7 + tmp8;
2562                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35;
2563                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0 + tmp1 + tmp11 + tmp12 + tmp13 + tmp14;
2564                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp9;
2565                                    }
2566                                }
2567                            } else { // constant data
2568                                for (index_t k=0; k<numEq; k++) {
2569                                    for (index_t m=0; m<numComp; m++) {
2570                                        const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
2571                                        const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
2572                                        const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
2573                                        const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
2574                                        const double tmp0 = 6*w1*(A_01 - A_10);
2575                                        const double tmp1 = 6*w1*(A_01 + A_10);
2576                                        const double tmp2 = 6*w1*(-A_01 - A_10);
2577                                        const double tmp3 = 6*w1*(-A_01 + A_10);
2578                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
2579                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
2580                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
2581                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
2582                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
2583                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
2584                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
2585                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
2586                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
2587                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
2588                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
2589                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
2590                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
2591                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
2592                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
2593                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
2594                                    }
2595                                }
2596                            }
2597                        }
2598                        ///////////////
2599                        // process B //
2600                        ///////////////
2601                        if (!B.isEmpty()) {
2602                            add_EM_S=true;
2603                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2604                            if (B.actsExpanded()) {
2605                                for (index_t k=0; k<numEq; k++) {
2606                                    for (index_t m=0; m<numComp; m++) {
2607                                        const double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
2608                                        const double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
2609                                        const double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
2610                                        const double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
2611                                        const double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];
2612                                        const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
2613                                        const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
2614                                        const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
2615                                        const double tmp0 = w11*(B_1_0 + B_1_1);
2616                                        const double tmp1 = w14*(B_1_2 + B_1_3);
2617                                        const double tmp2 = w15*(-B_0_1 - B_0_3);
2618                                        const double tmp3 = w10*(-B_0_0 - B_0_2);
2619                                        const double tmp4 = w11*(B_1_2 + B_1_3);
2620                                        const double tmp5 = w14*(B_1_0 + B_1_1);
2621                                        const double tmp6 = w11*(-B_1_2 - B_1_3);
2622                                        const double tmp7 = w14*(-B_1_0 - B_1_1);
2623                                        const double tmp8 = w11*(-B_1_0 - B_1_1);
2624                                        const double tmp9 = w14*(-B_1_2 - B_1_3);
2625                                        const double tmp10 = w10*(-B_0_1 - B_0_3);
2626                                        const double tmp11 = w15*(-B_0_0 - B_0_2);
2627                                        const double tmp12 = w15*(B_0_0 + B_0_2);
2628                                        const double tmp13 = w10*(B_0_1 + B_0_3);
2629                                        const double tmp14 = w10*(B_0_0 + B_0_2);
2630                                        const double tmp15 = w15*(B_0_1 + B_0_3);
2631                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=B_0_0*w12 + B_0_1*w10 + B_0_2*w15 + B_0_3*w13 + B_1_0*w16 + B_1_1*w14 + B_1_2*w11 + B_1_3*w17;
2632                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=B_0_0*w10 + B_0_1*w12 + B_0_2*w13 + B_0_3*w15 + tmp0 + tmp1;
2633                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=B_1_0*w11 + B_1_1*w17 + B_1_2*w16 + B_1_3*w14 + tmp14 + tmp15;
2634                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp4 + tmp5;
2635                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-B_0_0*w12 - B_0_1*w10 - B_0_2*w15 - B_0_3*w13 + tmp0 + tmp1;
2636                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-B_0_0*w10 - B_0_1*w12 - B_0_2*w13 - B_0_3*w15 + B_1_0*w14 + B_1_1*w16 + B_1_2*w17 + B_1_3*w11;
2637                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2638                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=B_1_0*w17 + B_1_1*w11 + B_1_2*w14 + B_1_3*w16 + tmp10 + tmp11;
2639                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-B_1_0*w16 - B_1_1*w14 - B_1_2*w11 - B_1_3*w17 + tmp14 + tmp15;
2640                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2641                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=B_0_0*w15 + B_0_1*w13 + B_0_2*w12 + B_0_3*w10 - B_1_0*w11 - B_1_1*w17 - B_1_2*w16 - B_1_3*w14;
2642                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=B_0_0*w13 + B_0_1*w15 + B_0_2*w10 + B_0_3*w12 + tmp6 + tmp7;
2643                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2 + tmp3 + tmp8 + tmp9;
2644                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-B_1_0*w14 - B_1_1*w16 - B_1_2*w17 - B_1_3*w11 + tmp10 + tmp11;
2645                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-B_0_0*w15 - B_0_1*w13 - B_0_2*w12 - B_0_3*w10 + tmp6 + tmp7;
2646                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-B_0_0*w13 - B_0_1*w15 - B_0_2*w10 - B_0_3*w12 - B_1_0*w17 - B_1_1*w11 - B_1_2*w14 - B_1_3*w16;
2647                                    }
2648                                }
2649                            } else { // constant data
2650                                for (index_t k=0; k<numEq; k++) {
2651                                    for (index_t m=0; m<numComp; m++) {
2652                                        const double wB0 = B_p[INDEX3(k,0,m,numEq,2)]*w18;
2653                                        const double wB1 = B_p[INDEX3(k,1,m,numEq,2)]*w19;
2654                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+= 2*wB0 + 2*wB1;
2655                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= 2*wB0 +   wB1;
2656                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=   wB0 + 2*wB1;
2657                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=   wB0 +   wB1;
2658                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-2*wB0 +   wB1;
2659                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-2*wB0 + 2*wB1;
2660                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=  -wB0 +   wB1;
2661                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=  -wB0 + 2*wB1;
2662                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=   wB0 - 2*wB1;
2663                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=   wB0 -   wB1;
2664                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= 2*wB0 - 2*wB1;
2665                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= 2*wB0 -   wB1;
2666                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=  -wB0 -   wB1;
2667                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=  -wB0 - 2*wB1;
2668                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-2*wB0 -   wB1;
2669                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-2*wB0 - 2*wB1;
2670                                    }
2671                                }
2672                            }
2673                        }
2674                        ///////////////
2675                        // process C //
2676                        ///////////////
2677                        if (!C.isEmpty()) {
2678                            add_EM_S=true;
2679                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2680                            if (C.actsExpanded()) {
2681                                for (index_t k=0; k<numEq; k++) {
2682                                    for (index_t m=0; m<numComp; m++) {
2683                                        const double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
2684                                        const double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
2685                                        const double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
2686                                        const double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
2687                                        const double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];
2688                                        const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
2689                               &nbs