/[escript]/branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp

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

revision 3722 by caltinay, Wed Dec 7 05:53:22 2011 UTC revision 3766 by caltinay, Thu Jan 12 08:17:49 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2011 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 13  Line 13 
13    
14  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
15  extern "C" {  extern "C" {
16  #include "paso/SystemMatrixPattern.h"  #include "paso/SystemMatrix.h"
17  }  }
18    
19  #if USE_SILO  #if USE_SILO
# Line 43  Rectangle::Rectangle(int n0, int n1, dou Line 43  Rectangle::Rectangle(int n0, int n1, dou
43      if (m_NX*m_NY != m_mpiInfo->size)      if (m_NX*m_NY != m_mpiInfo->size)
44          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
45    
46      if (n0%m_NX > 0 || n1%m_NY > 0)      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)
47          throw RipleyException("Number of elements must be separable into number of ranks in each dimension");          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");
48    
49      // local number of elements      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
50      m_NE0 = n0/m_NX;          throw RipleyException("Too few elements for the number of ranks");
51      m_NE1 = n1/m_NY;  
52      // local number of nodes (not necessarily owned)      // local number of elements (with and without overlap)
53        m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
54        if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
55            m_NE0++;
56        else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)
57            m_ownNE0--;
58    
59        m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
60        if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)
61            m_NE1++;
62        else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)
63            m_ownNE1--;
64    
65        // local number of nodes
66      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
67      m_N1 = m_NE1+1;      m_N1 = m_NE1+1;
68    
69      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
70      m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);
71      m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);      if (m_offset0 > 0)
72            m_offset0--;
73        m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);
74        if (m_offset1 > 0)
75            m_offset1--;
76    
77      populateSampleIds();      populateSampleIds();
78        createPattern();
79    
80        if (m_gNE0>2*m_NX && m_gNE1>2*m_NY && m_gNE0%2==0 && m_gNE1%2==0) {
81            m_coarseMesh.reset(new Rectangle(n0/2, n1/2, l0, l1, d0, d1));
82        }
83  }  }
84    
85  Rectangle::~Rectangle()  Rectangle::~Rectangle()
# Line 69  string Rectangle::getDescription() const Line 93  string Rectangle::getDescription() const
93    
94  bool Rectangle::operator==(const AbstractDomain& other) const  bool Rectangle::operator==(const AbstractDomain& other) const
95  {  {
96      if (dynamic_cast<const Rectangle*>(&other))      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
97          return this==&other;      if (o) {
98            return (RipleyDomain::operator==(other) &&
99                    m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1
100                    && m_l0==o->m_l0 && m_l1==o->m_l1
101                    && m_NX==o->m_NX && m_NY==o->m_NY);
102        }
103    
104      return false;      return false;
105  }  }
# Line 86  void Rectangle::dump(const string& fileN Line 115  void Rectangle::dump(const string& fileN
115      const int NUM_SILO_FILES = 1;      const int NUM_SILO_FILES = 1;
116      const char* blockDirFmt = "/block%04d";      const char* blockDirFmt = "/block%04d";
117      int driver=DB_HDF5;          int driver=DB_HDF5;    
     string siloPath;  
118      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
119    
120  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 106  void Rectangle::dump(const string& fileN Line 134  void Rectangle::dump(const string& fileN
134                          PMPIO_DefaultClose, (void*)&driver);                          PMPIO_DefaultClose, (void*)&driver);
135          }          }
136          if (baton) {          if (baton) {
137              char str[64];              char siloPath[64];
138              snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));              snprintf(siloPath, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
139              siloPath = str;              dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath);
             dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());  
140          }          }
141  #endif  #endif
142      } else {      } else {
# Line 121  void Rectangle::dump(const string& fileN Line 148  void Rectangle::dump(const string& fileN
148              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
149                      getDescription().c_str(), driver);                      getDescription().c_str(), driver);
150          }          }
151            char siloPath[64];
152            snprintf(siloPath, 64, blockDirFmt, 0);
153            DBMkDir(dbfile, siloPath);
154            DBSetDir(dbfile, siloPath);
155      }      }
156    
157      if (!dbfile)      if (!dbfile)
# Line 222  const int* Rectangle::borrowSampleRefere Line 253  const int* Rectangle::borrowSampleRefere
253      switch (fsType) {      switch (fsType) {
254          case Nodes:          case Nodes:
255              return &m_nodeId[0];              return &m_nodeId[0];
256            case ReducedNodes:
257                if (m_coarseMesh)
258                    return m_coarseMesh->borrowSampleReferenceIDs(Nodes);
259                break;
260            case DegreesOfFreedom:
261                return &m_dofId[0];
262            case ReducedDegreesOfFreedom:
263                if (m_coarseMesh)
264                    return m_coarseMesh->borrowSampleReferenceIDs(DegreesOfFreedom);
265                break;
266          case Elements:          case Elements:
267            case ReducedElements:
268              return &m_elementId[0];              return &m_elementId[0];
269          case FaceElements:          case FaceElements:
270            case ReducedFaceElements:
271              return &m_faceId[0];              return &m_faceId[0];
272          default:          default:
273              break;              break;
# Line 236  const int* Rectangle::borrowSampleRefere Line 279  const int* Rectangle::borrowSampleRefere
279      throw RipleyException(msg.str());      throw RipleyException(msg.str());
280  }  }
281    
282  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsType, index_t id) const
283  {  {
284  #ifdef ESYS_MPI      if (getMPISize()==1)
285      if (fsCode != ReducedNodes) {          return true;
286          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];  
287          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;      switch (fsType) {
288          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);          case Nodes:
289                return (m_dofMap[id] < getNumDOF());
290            case ReducedNodes:
291                if (m_coarseMesh)
292                    return m_coarseMesh->ownSample(Nodes, id);
293                break;
294            case DegreesOfFreedom:
295            case ReducedDegreesOfFreedom:
296                return true;
297            case Elements:
298            case ReducedElements:
299                // check ownership of element's bottom left node
300                return (m_dofMap[id%m_NE0+m_N0*(id/m_NE0)] < getNumDOF());
301            case FaceElements:
302            case ReducedFaceElements:
303                {
304                    // determine which face the sample belongs to before
305                    // checking ownership of face element's first node
306                    const IndexVector faces = getNumFacesPerBoundary();
307                    dim_t n=0;
308                    for (size_t i=0; i<faces.size(); i++) {
309                        n+=faces[i];
310                        if (id<n) {
311                            index_t k;
312                            if (i==1)
313                                k=m_N0-1;
314                            else if (i==3)
315                                k=m_N0*(m_N1-1);
316                            else
317                                k=0;
318                            // determine whether to move right or up
319                            const index_t delta=(i/2==0 ? m_N0 : 1);
320                            return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());
321                        }
322                    }
323                    return false;
324                }
325            default:
326                break;
327        }
328    
329        stringstream msg;
330        msg << "ownSample() not implemented for "
331            << functionSpaceTypeAsString(fsType);
332        throw RipleyException(msg.str());
333    }
334    
335    void Rectangle::setToNormal(escript::Data& out) const
336    {
337        if (out.getFunctionSpace().getTypeCode() == FaceElements) {
338            out.requireWrite();
339    #pragma omp parallel
340            {
341                if (m_faceOffset[0] > -1) {
342    #pragma omp for nowait
343                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
344                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
345                        // set vector at two quadrature points
346                        *o++ = -1.;
347                        *o++ = 0.;
348                        *o++ = -1.;
349                        *o = 0.;
350                    }
351                }
352    
353                if (m_faceOffset[1] > -1) {
354    #pragma omp for nowait
355                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
356                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
357                        // set vector at two quadrature points
358                        *o++ = 1.;
359                        *o++ = 0.;
360                        *o++ = 1.;
361                        *o = 0.;
362                    }
363                }
364    
365                if (m_faceOffset[2] > -1) {
366    #pragma omp for nowait
367                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
368                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
369                        // set vector at two quadrature points
370                        *o++ = 0.;
371                        *o++ = -1.;
372                        *o++ = 0.;
373                        *o = -1.;
374                    }
375                }
376    
377                if (m_faceOffset[3] > -1) {
378    #pragma omp for nowait
379                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
380                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
381                        // set vector at two quadrature points
382                        *o++ = 0.;
383                        *o++ = 1.;
384                        *o++ = 0.;
385                        *o = 1.;
386                    }
387                }
388            } // end of parallel section
389        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
390            out.requireWrite();
391    #pragma omp parallel
392            {
393                if (m_faceOffset[0] > -1) {
394    #pragma omp for nowait
395                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
396                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
397                        *o++ = -1.;
398                        *o = 0.;
399                    }
400                }
401    
402                if (m_faceOffset[1] > -1) {
403    #pragma omp for nowait
404                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
405                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
406                        *o++ = 1.;
407                        *o = 0.;
408                    }
409                }
410    
411                if (m_faceOffset[2] > -1) {
412    #pragma omp for nowait
413                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
414                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
415                        *o++ = 0.;
416                        *o = -1.;
417                    }
418                }
419    
420                if (m_faceOffset[3] > -1) {
421    #pragma omp for nowait
422                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
423                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
424                        *o++ = 0.;
425                        *o = 1.;
426                    }
427                }
428            } // end of parallel section
429    
430      } else {      } else {
431          stringstream msg;          stringstream msg;
432          msg << "ownSample() not implemented for "          msg << "setToNormal() not implemented for "
433              << functionSpaceTypeAsString(fsCode);              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
434          throw RipleyException(msg.str());          throw RipleyException(msg.str());
435      }      }
 #else  
     return true;  
 #endif  
436  }  }
437    
438  void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const  void Rectangle::setToSize(escript::Data& out) const
439    {
440        if (out.getFunctionSpace().getTypeCode() == Elements
441                || out.getFunctionSpace().getTypeCode() == ReducedElements) {
442            out.requireWrite();
443            const dim_t numQuad=out.getNumDataPointsPerSample();
444            const double hSize=getFirstCoordAndSpacing(0).second;
445            const double vSize=getFirstCoordAndSpacing(1).second;
446            const double size=min(hSize,vSize);
447    #pragma omp parallel for
448            for (index_t k = 0; k < getNumElements(); ++k) {
449                double* o = out.getSampleDataRW(k);
450                fill(o, o+numQuad, size);
451            }
452        } else if (out.getFunctionSpace().getTypeCode() == FaceElements
453                || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
454            out.requireWrite();
455            const dim_t numQuad=out.getNumDataPointsPerSample();
456            const double hSize=getFirstCoordAndSpacing(0).second;
457            const double vSize=getFirstCoordAndSpacing(1).second;
458    #pragma omp parallel
459            {
460                if (m_faceOffset[0] > -1) {
461    #pragma omp for nowait
462                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
463                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
464                        fill(o, o+numQuad, vSize);
465                    }
466                }
467    
468                if (m_faceOffset[1] > -1) {
469    #pragma omp for nowait
470                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
471                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
472                        fill(o, o+numQuad, vSize);
473                    }
474                }
475    
476                if (m_faceOffset[2] > -1) {
477    #pragma omp for nowait
478                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
479                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
480                        fill(o, o+numQuad, hSize);
481                    }
482                }
483    
484                if (m_faceOffset[3] > -1) {
485    #pragma omp for nowait
486                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
487                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
488                        fill(o, o+numQuad, hSize);
489                    }
490                }
491            } // end of parallel section
492    
493        } else {
494            stringstream msg;
495            msg << "setToSize() not implemented for "
496                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
497            throw RipleyException(msg.str());
498        }
499    }
500    
501    Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
502                                                    bool reducedColOrder) const
503    {
504        if (reducedRowOrder || reducedColOrder)
505            throw RipleyException("getPattern() not implemented for reduced order");
506    
507        return m_pattern;
508    }
509    
510    void Rectangle::Print_Mesh_Info(const bool full) const
511    {
512        RipleyDomain::Print_Mesh_Info(full);
513        if (full) {
514            cout << "     Id  Coordinates" << endl;
515            cout.precision(15);
516            cout.setf(ios::scientific, ios::floatfield);
517            pair<double,double> xdx = getFirstCoordAndSpacing(0);
518            pair<double,double> ydy = getFirstCoordAndSpacing(1);
519            for (index_t i=0; i < getNumNodes(); i++) {
520                cout << "  " << setw(5) << m_nodeId[i]
521                    << "  " << xdx.first+(i%m_N0)*xdx.second
522                    << "  " << ydy.first+(i/m_N0)*ydy.second << endl;
523            }
524        }
525    }
526    
527    IndexVector Rectangle::getNumNodesPerDim() const
528    {
529        IndexVector ret;
530        ret.push_back(m_N0);
531        ret.push_back(m_N1);
532        return ret;
533    }
534    
535    IndexVector Rectangle::getNumElementsPerDim() const
536    {
537        IndexVector ret;
538        ret.push_back(m_NE0);
539        ret.push_back(m_NE1);
540        return ret;
541    }
542    
543    IndexVector Rectangle::getNumFacesPerBoundary() const
544    {
545        IndexVector ret(4, 0);
546        //left
547        if (m_offset0==0)
548            ret[0]=m_NE1;
549        //right
550        if (m_mpiInfo->rank%m_NX==m_NX-1)
551            ret[1]=m_NE1;
552        //bottom
553        if (m_offset1==0)
554            ret[2]=m_NE0;
555        //top
556        if (m_mpiInfo->rank/m_NX==m_NY-1)
557            ret[3]=m_NE0;
558        return ret;
559    }
560    
561    IndexVector Rectangle::getNumSubdivisionsPerDim() const
562    {
563        IndexVector ret;
564        ret.push_back(m_NX);
565        ret.push_back(m_NY);
566        return ret;
567    }
568    
569    pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
570    {
571        if (dim==0) {
572            return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
573        } else if (dim==1) {
574            return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
575        }
576        throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
577    }
578    
579    //protected
580    dim_t Rectangle::getNumDOF() const
581    {
582        return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;
583    }
584    
585    //protected
586    dim_t Rectangle::getNumFaceElements() const
587    {
588        const IndexVector faces = getNumFacesPerBoundary();
589        dim_t n=0;
590        for (size_t i=0; i<faces.size(); i++)
591            n+=faces[i];
592        return n;
593    }
594    
595    //protected
596    void Rectangle::assembleCoordinates(escript::Data& arg) const
597    {
598        escriptDataC x = arg.getDataC();
599        int numDim = m_numDim;
600        if (!isDataPointShapeEqual(&x, 1, &numDim))
601            throw RipleyException("setToX: Invalid Data object shape");
602        if (!numSamplesEqual(&x, 1, getNumNodes()))
603            throw RipleyException("setToX: Illegal number of samples in Data object");
604    
605        pair<double,double> xdx = getFirstCoordAndSpacing(0);
606        pair<double,double> ydy = getFirstCoordAndSpacing(1);
607        arg.requireWrite();
608    #pragma omp parallel for
609        for (dim_t i1 = 0; i1 < m_N1; i1++) {
610            for (dim_t i0 = 0; i0 < m_N0; i0++) {
611                double* point = arg.getSampleDataRW(i0+m_N0*i1);
612                point[0] = xdx.first+i0*xdx.second;
613                point[1] = ydy.first+i1*ydy.second;
614            }
615        }
616    }
617    
618    //protected
619    void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
620  {  {
     escript::Data& in = *const_cast<escript::Data*>(&cIn);  
621      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
622      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
623      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
624        const double cx0 = -1./h0;
625        const double cx1 = -.78867513459481288225/h0;
626        const double cx2 = -.5/h0;
627        const double cx3 = -.21132486540518711775/h0;
628        const double cx4 = .21132486540518711775/h0;
629        const double cx5 = .5/h0;
630        const double cx6 = .78867513459481288225/h0;
631        const double cx7 = 1./h0;
632        const double cy0 = -1./h1;
633        const double cy1 = -.78867513459481288225/h1;
634        const double cy2 = -.5/h1;
635        const double cy3 = -.21132486540518711775/h1;
636        const double cy4 = .21132486540518711775/h1;
637        const double cy5 = .5/h1;
638        const double cy6 = .78867513459481288225/h1;
639        const double cy7 = 1./h1;
640    
641      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
642          /* GENERATOR SNIP_GRAD_ELEMENTS TOP */          out.requireWrite();
         const double tmp0_13 = 0.78867513459481288225/h1;  
         const double tmp0_0 = 0.78867513459481288225/h0;  
         const double tmp0_4 = 0.21132486540518711775/h0;  
         const double tmp0_10 = 0.78867513459481288225/h1;  
         const double tmp0_1 = 0.21132486540518711775/h0;  
         const double tmp0_8 = -0.21132486540518711775/h1;  
         const double tmp0_14 = 0.21132486540518711775/h1;  
         const double tmp0_5 = 0.78867513459481288225/h0;  
         const double tmp0_11 = -0.78867513459481288225/h1;  
         const double tmp0_2 = -0.21132486540518711775/h0;  
         const double tmp0_9 = 0.21132486540518711775/h1;  
         const double tmp0_15 = -0.21132486540518711775/h1;  
         const double tmp0_6 = -0.78867513459481288225/h0;  
         const double tmp0_3 = -0.78867513459481288225/h0;  
         const double tmp0_12 = -0.78867513459481288225/h1;  
         const double tmp0_7 = -0.21132486540518711775/h0;  
643  #pragma omp parallel for  #pragma omp parallel for
644          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
645              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 287  void Rectangle::setToGradient(escript::D Line 649  void Rectangle::setToGradient(escript::D
649                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
650                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
651                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
652                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2 + f_10[i]*tmp0_0 + f_11[i]*tmp0_1;                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
653                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_11 + f_01[i]*tmp0_10 + f_10[i]*tmp0_8 + f_11[i]*tmp0_9;                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
654                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2 + f_10[i]*tmp0_0 + f_11[i]*tmp0_1;                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
655                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_15 + f_01[i]*tmp0_14 + f_10[i]*tmp0_12 + f_11[i]*tmp0_13;                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
656                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
657                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*tmp0_11 + f_01[i]*tmp0_10 + f_10[i]*tmp0_8 + f_11[i]*tmp0_9;                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
658                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
659                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*tmp0_15 + f_01[i]*tmp0_14 + f_10[i]*tmp0_12 + f_11[i]*tmp0_13;                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
660                  } /* end of component loop i */                  } /* end of component loop i */
661              } /* end of k0 loop */              } /* end of k0 loop */
662          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */  
663      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
664          /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */          out.requireWrite();
         const double tmp0_3 = 0.5/h1;  
         const double tmp0_2 = -0.5/h1;  
         const double tmp0_1 = -0.5/h0;  
         const double tmp0_0 = 0.5/h0;  
665  #pragma omp parallel for  #pragma omp parallel for
666          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
667              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 314  void Rectangle::setToGradient(escript::D Line 671  void Rectangle::setToGradient(escript::D
671                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
672                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
673                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
674                      o[INDEX3(i,0,0,numComp,2)] = tmp0_0*(f_10[i] + f_11[i]) + tmp0_1*(f_00[i] + f_01[i]);                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
675                      o[INDEX3(i,1,0,numComp,2)] = tmp0_2*(f_00[i] + f_10[i]) + tmp0_3*(f_01[i] + f_11[i]);                      o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
676                  } /* end of component loop i */                  } /* end of component loop i */
677              } /* end of k0 loop */              } /* end of k0 loop */
678          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */  
679      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
680            out.requireWrite();
681  #pragma omp parallel  #pragma omp parallel
682          {          {
             /* GENERATOR SNIP_GRAD_FACES TOP */  
683              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
                 const double tmp0_0 = 0.78867513459481288225/h0;  
                 const double tmp0_4 = 0.21132486540518711775/h0;  
                 const double tmp0_1 = 0.21132486540518711775/h0;  
                 const double tmp0_8 = 1.0/h1;  
                 const double tmp0_5 = 0.78867513459481288225/h0;  
                 const double tmp0_2 = -0.21132486540518711775/h0;  
                 const double tmp0_9 = -1/h1;  
                 const double tmp0_6 = -0.78867513459481288225/h0;  
                 const double tmp0_3 = -0.78867513459481288225/h0;  
                 const double tmp0_7 = -0.21132486540518711775/h0;  
684  #pragma omp for nowait  #pragma omp for nowait
685                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
686                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
# Line 343  void Rectangle::setToGradient(escript::D Line 689  void Rectangle::setToGradient(escript::D
689                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
690                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
691                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
692                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2 + f_10[i]*tmp0_0 + f_11[i]*tmp0_1;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
693                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
694                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
695                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
696                      } /* end of component loop i */                      } /* end of component loop i */
697                  } /* end of k1 loop */                  } /* end of k1 loop */
698              } /* end of face 0 */              } /* end of face 0 */
699              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
                 const double tmp0_0 = 0.78867513459481288225/h0;  
                 const double tmp0_4 = 0.21132486540518711775/h0;  
                 const double tmp0_1 = 0.21132486540518711775/h0;  
                 const double tmp0_8 = -1/h1;  
                 const double tmp0_5 = 0.78867513459481288225/h0;  
                 const double tmp0_2 = -0.21132486540518711775/h0;  
                 const double tmp0_9 = 1.0/h1;  
                 const double tmp0_6 = -0.78867513459481288225/h0;  
                 const double tmp0_3 = -0.78867513459481288225/h0;  
                 const double tmp0_7 = -0.21132486540518711775/h0;  
700  #pragma omp for nowait  #pragma omp for nowait
701                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
702                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
# Line 369  void Rectangle::setToGradient(escript::D Line 705  void Rectangle::setToGradient(escript::D
705                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
706                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
707                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
708                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2 + f_10[i]*tmp0_0 + f_11[i]*tmp0_1;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
709                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
710                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
711                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
712                      } /* end of component loop i */                      } /* end of component loop i */
713                  } /* end of k1 loop */                  } /* end of k1 loop */
714              } /* end of face 1 */              } /* end of face 1 */
715              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
                 const double tmp0_0 = -1/h0;  
                 const double tmp0_4 = 0.21132486540518711775/h1;  
                 const double tmp0_1 = 1.0/h0;  
                 const double tmp0_8 = 0.78867513459481288225/h1;  
                 const double tmp0_5 = 0.78867513459481288225/h1;  
                 const double tmp0_2 = -0.78867513459481288225/h1;  
                 const double tmp0_9 = 0.21132486540518711775/h1;  
                 const double tmp0_6 = -0.21132486540518711775/h1;  
                 const double tmp0_3 = -0.21132486540518711775/h1;  
                 const double tmp0_7 = -0.78867513459481288225/h1;  
716  #pragma omp for nowait  #pragma omp for nowait
717                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
718                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
# Line 395  void Rectangle::setToGradient(escript::D Line 721  void Rectangle::setToGradient(escript::D
721                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
722                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
723                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
724                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
725                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_2 + f_01[i]*tmp0_5 + f_10[i]*tmp0_3 + f_11[i]*tmp0_4;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
726                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
727                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_6 + f_01[i]*tmp0_9 + f_10[i]*tmp0_7 + f_11[i]*tmp0_8;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
728                      } /* end of component loop i */                      } /* end of component loop i */
729                  } /* end of k0 loop */                  } /* end of k0 loop */
730              } /* end of face 2 */              } /* end of face 2 */
731              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
                 const double tmp0_0 = 1.0/h0;  
                 const double tmp0_4 = -0.21132486540518711775/h1;  
                 const double tmp0_1 = -1/h0;  
                 const double tmp0_8 = -0.78867513459481288225/h1;  
                 const double tmp0_5 = -0.78867513459481288225/h1;  
                 const double tmp0_2 = 0.21132486540518711775/h1;  
                 const double tmp0_9 = -0.21132486540518711775/h1;  
                 const double tmp0_6 = 0.78867513459481288225/h1;  
                 const double tmp0_3 = 0.78867513459481288225/h1;  
                 const double tmp0_7 = 0.21132486540518711775/h1;  
732  #pragma omp for nowait  #pragma omp for nowait
733                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
734                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
# Line 421  void Rectangle::setToGradient(escript::D Line 737  void Rectangle::setToGradient(escript::D
737                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
738                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
739                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
740                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
741                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_5 + f_01[i]*tmp0_3 + f_10[i]*tmp0_4 + f_11[i]*tmp0_2;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
742                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
743                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_7 + f_10[i]*tmp0_8 + f_11[i]*tmp0_6;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
744                      } /* end of component loop i */                      } /* end of component loop i */
745                  } /* end of k0 loop */                  } /* end of k0 loop */
746              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_GRAD_FACES BOTTOM */  
747          } // end of parallel section          } // end of parallel section
748      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
749            out.requireWrite();
750  #pragma omp parallel  #pragma omp parallel
751          {          {
             /* GENERATOR SNIP_GRAD_REDUCED_FACES TOP */  
752              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
                 const double tmp0_3 = -1/h1;  
                 const double tmp0_2 = 1.0/h1;  
                 const double tmp0_1 = -0.5/h0;  
                 const double tmp0_0 = 0.5/h0;  
753  #pragma omp for nowait  #pragma omp for nowait
754                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
755                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
# Line 447  void Rectangle::setToGradient(escript::D Line 758  void Rectangle::setToGradient(escript::D
758                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
759                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
760                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
761                          o[INDEX3(i,0,0,numComp,2)] = tmp0_0*(f_10[i] + f_11[i]) + tmp0_1*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
762                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
763                      } /* end of component loop i */                      } /* end of component loop i */
764                  } /* end of k1 loop */                  } /* end of k1 loop */
765              } /* end of face 0 */              } /* end of face 0 */
766              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
                 const double tmp0_3 = 1.0/h1;  
                 const double tmp0_2 = -1/h1;  
                 const double tmp0_1 = -0.5/h0;  
                 const double tmp0_0 = 0.5/h0;  
767  #pragma omp for nowait  #pragma omp for nowait
768                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
769                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
# Line 465  void Rectangle::setToGradient(escript::D Line 772  void Rectangle::setToGradient(escript::D
772                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
773                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
774                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
775                          o[INDEX3(i,0,0,numComp,2)] = tmp0_0*(f_10[i] + f_11[i]) + tmp0_1*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
776                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*tmp0_2 + f_11[i]*tmp0_3;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
777                      } /* end of component loop i */                      } /* end of component loop i */
778                  } /* end of k1 loop */                  } /* end of k1 loop */
779              } /* end of face 1 */              } /* end of face 1 */
780              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
                 const double tmp0_3 = 0.5/h1;  
                 const double tmp0_2 = -0.5/h1;  
                 const double tmp0_1 = 1.0/h0;  
                 const double tmp0_0 = -1/h0;  
781  #pragma omp for nowait  #pragma omp for nowait
782                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
783                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
# Line 483  void Rectangle::setToGradient(escript::D Line 786  void Rectangle::setToGradient(escript::D
786                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
787                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
788                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
789                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
790                          o[INDEX3(i,1,0,numComp,2)] = tmp0_2*(f_00[i] + f_10[i]) + tmp0_3*(f_01[i] + f_11[i]);                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
791                      } /* end of component loop i */                      } /* end of component loop i */
792                  } /* end of k0 loop */                  } /* end of k0 loop */
793              } /* end of face 2 */              } /* end of face 2 */
794              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
                 const double tmp0_3 = -0.5/h1;  
                 const double tmp0_2 = 0.5/h1;  
                 const double tmp0_1 = -1/h0;  
                 const double tmp0_0 = 1.0/h0;  
795  #pragma omp for nowait  #pragma omp for nowait
796                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
797                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
# Line 501  void Rectangle::setToGradient(escript::D Line 800  void Rectangle::setToGradient(escript::D
800                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
801                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
802                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
803                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
804                          o[INDEX3(i,1,0,numComp,2)] = tmp0_2*(f_01[i] + f_11[i]) + tmp0_3*(f_00[i] + f_10[i]);                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
805                      } /* end of component loop i */                      } /* end of component loop i */
806                  } /* end of k0 loop */                  } /* end of k0 loop */
807              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */  
808          } // end of parallel section          } // end of parallel section
     } else {  
         stringstream msg;  
         msg << "setToGradient() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
809      }      }
810  }  }
811    
812  void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  //protected
813    void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
814  {  {
815      escript::Data& in = *const_cast<escript::Data*>(&arg);      const dim_t numComp = arg.getDataPointSize();
     const dim_t numComp = in.getDataPointSize();  
816      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
817      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
818        const index_t left = (m_offset0==0 ? 0 : 1);
819        const index_t bottom = (m_offset1==0 ? 0 : 1);
820      if (arg.getFunctionSpace().getTypeCode() == Elements) {      if (arg.getFunctionSpace().getTypeCode() == Elements) {
821          const double w_0 = h0*h1/4.;          const double w = h0*h1/4.;
822  #pragma omp parallel  #pragma omp parallel
823          {          {
824              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
825  #pragma omp for nowait  #pragma omp for nowait
826              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
827                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
828                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
829                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
830                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
831                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
832                          const register double f_2 = f[INDEX2(i,2,numComp)];                          const register double f2 = f[INDEX2(i,2,numComp)];
833                          const register double f_3 = f[INDEX2(i,3,numComp)];                          const register double f3 = f[INDEX2(i,3,numComp)];
834                          int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;                          int_local[i]+=(f0+f1+f2+f3)*w;
835                      }  /* end of component loop i */                      }  /* end of component loop i */
836                  } /* end of k0 loop */                  } /* end of k0 loop */
837              } /* end of k1 loop */              } /* end of k1 loop */
# Line 546  void Rectangle::setToIntegrals(vector<do Line 841  void Rectangle::setToIntegrals(vector<do
841                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
842          } // end of parallel section          } // end of parallel section
843      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
844          const double w_0 = h0*h1;          const double w = h0*h1;
845  #pragma omp parallel  #pragma omp parallel
846          {          {
847              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
848  #pragma omp for nowait  #pragma omp for nowait
849              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
850                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
851                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
852                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
853                          int_local[i]+=f[i]*w_0;                          int_local[i]+=f[i]*w;
854                      }  /* end of component loop i */                      }  /* end of component loop i */
855                  } /* end of k0 loop */                  } /* end of k0 loop */
856              } /* end of k1 loop */              } /* end of k1 loop */
# Line 565  void Rectangle::setToIntegrals(vector<do Line 860  void Rectangle::setToIntegrals(vector<do
860                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
861          } // end of parallel section          } // end of parallel section
862      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
863          const double w_0 = h0/2.;          const double w0 = h0/2.;
864          const double w_1 = h1/2.;          const double w1 = h1/2.;
865  #pragma omp parallel  #pragma omp parallel
866          {          {
867              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
868              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
869  #pragma omp for nowait  #pragma omp for nowait
870                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
871                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
872                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
873                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
874                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
875                          int_local[i]+=(f_0+f_1)*w_1;                          int_local[i]+=(f0+f1)*w1;
876                      }  /* end of component loop i */                      }  /* end of component loop i */
877                  } /* end of k1 loop */                  } /* end of k1 loop */
878              }              }
879    
880              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
881  #pragma omp for nowait  #pragma omp for nowait
882                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
883                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
884                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
885                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
886                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
887                          int_local[i]+=(f_0+f_1)*w_1;                          int_local[i]+=(f0+f1)*w1;
888                      }  /* end of component loop i */                      }  /* end of component loop i */
889                  } /* end of k1 loop */                  } /* end of k1 loop */
890              }              }
891    
892              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
893  #pragma omp for nowait  #pragma omp for nowait
894                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
895                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
896                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
897                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
898                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
899                          int_local[i]+=(f_0+f_1)*w_0;                          int_local[i]+=(f0+f1)*w0;
900                      }  /* end of component loop i */                      }  /* end of component loop i */
901                  } /* end of k0 loop */                  } /* end of k0 loop */
902              }              }
903    
904              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
905  #pragma omp for nowait  #pragma omp for nowait
906                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
907                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
908                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
909                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
910                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
911                          int_local[i]+=(f_0+f_1)*w_0;                          int_local[i]+=(f0+f1)*w0;
912                      }  /* end of component loop i */                      }  /* end of component loop i */
913                  } /* end of k0 loop */                  } /* end of k0 loop */
914              }              }
# Line 628  void Rectangle::setToIntegrals(vector<do Line 923  void Rectangle::setToIntegrals(vector<do
923              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
924              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
925  #pragma omp for nowait  #pragma omp for nowait
926                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
927                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
928                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
929                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
930                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 638  void Rectangle::setToIntegrals(vector<do Line 933  void Rectangle::setToIntegrals(vector<do
933    
934              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
935  #pragma omp for nowait  #pragma omp for nowait
936                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
937                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
938                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
939                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
940                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 648  void Rectangle::setToIntegrals(vector<do Line 943  void Rectangle::setToIntegrals(vector<do
943    
944              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
945  #pragma omp for nowait  #pragma omp for nowait
946                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
947                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
948                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
949                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
950                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 658  void Rectangle::setToIntegrals(vector<do Line 953  void Rectangle::setToIntegrals(vector<do
953    
954              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
955  #pragma omp for nowait  #pragma omp for nowait
956                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
957                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
958                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
959                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
960                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 670  void Rectangle::setToIntegrals(vector<do Line 965  void Rectangle::setToIntegrals(vector<do
965              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
966                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
967          } // end of parallel section          } // end of parallel section
     } else {  
         stringstream msg;  
         msg << "setToIntegrals() not implemented for "  
             << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
968      }      }
969  }  }
970    
971  void Rectangle::setToNormal(escript::Data& out) const  //protected
972    dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
973  {  {
974      if (out.getFunctionSpace().getTypeCode() == FaceElements) {      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
975  #pragma omp parallel      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
976          {      const int x=node%nDOF0;
977              if (m_faceOffset[0] > -1) {      const int y=node/nDOF0;
978  #pragma omp for nowait      dim_t num=0;
979                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {      // loop through potential neighbours and add to index if positions are
980                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);      // within bounds
981                      // set vector at two quadrature points      for (int i1=-1; i1<2; i1++) {
982                      *o++ = -1.;          for (int i0=-1; i0<2; i0++) {
983                      *o++ = 0.;              // skip node itself
984                      *o++ = -1.;              if (i0==0 && i1==0)
985                      *o = 0.;                  continue;
986                  }              // location of neighbour node
987              }              const int nx=x+i0;
988                const int ny=y+i1;
989              if (m_faceOffset[1] > -1) {              if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {
990  #pragma omp for nowait                  index.push_back(ny*nDOF0+nx);
991                  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.;  
                 }  
             }  
   
             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.;  
                 }  
992              }              }
993            }
             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());  
994      }      }
 }  
995    
996  void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,      return num;
         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");  
997  }  }
998    
999  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  //protected
1000                                                  bool reducedColOrder) const  void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const
1001  {  {
1002      if (reducedRowOrder || reducedColOrder)      const dim_t numComp = in.getDataPointSize();
1003          throw RipleyException("getPattern() not implemented for reduced order");      out.requireWrite();
1004    
     // connector  
     RankVector neighbour;  
     IndexVector offsetInShared(1,0);  
     IndexVector sendShared, recvShared;  
     const IndexVector faces=getNumFacesPerBoundary();  
1005      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1006      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1007      // corner node from bottom-left      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1008      if (faces[0] == 0 && faces[2] == 0) {      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1009          neighbour.push_back(m_mpiInfo->rank-m_NX-1);  #pragma omp parallel for
1010          offsetInShared.push_back(offsetInShared.back()+1);      for (index_t i=0; i<nDOF1; i++) {
1011          sendShared.push_back(m_nodeId[m_N0+left]);          for (index_t j=0; j<nDOF0; j++) {
1012          recvShared.push_back(m_nodeId[0]);              const index_t n=j+left+(i+bottom)*m_N0;
1013      }              const double* src=in.getSampleDataRO(n);
1014      // bottom edge              copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
     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;  
1015          }          }
1016      }      }
1017  }  }
1018    
 IndexVector Rectangle::getNumNodesPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_N0);  
     ret.push_back(m_N1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumElementsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NE0);  
     ret.push_back(m_NE1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumFacesPerBoundary() const  
 {  
     IndexVector ret(4, 0);  
     //left  
     if (m_offset0==0)  
         ret[0]=m_NE1;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         ret[1]=m_NE1;  
     //bottom  
     if (m_offset1==0)  
         ret[2]=m_NE0;  
     //top  
     if (m_mpiInfo->rank/m_NX==m_NY-1)  
         ret[3]=m_NE0;  
     return ret;  
 }  
   
 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");  
 }  
   
 //protected  
 dim_t Rectangle::getNumFaceElements() const  
 {  
     const IndexVector faces = getNumFacesPerBoundary();  
     dim_t n=0;  
     for (size_t i=0; i<faces.size(); i++)  
         n+=faces[i];  
     return n;  
 }  
   
1019  //protected  //protected
1020  void Rectangle::assembleCoordinates(escript::Data& arg) const  void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const
1021  {  {
1022      escriptDataC x = arg.getDataC();      const dim_t numComp = in.getDataPointSize();
1023      int numDim = m_numDim;      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
1024      if (!isDataPointShapeEqual(&x, 1, &numDim))      in.requireWrite();
1025          throw RipleyException("setToX: Invalid Data object shape");      Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));
1026      if (!numSamplesEqual(&x, 1, getNumNodes()))  
1027          throw RipleyException("setToX: Illegal number of samples in Data object");      const dim_t numDOF = getNumDOF();
1028        out.requireWrite();
1029        const double* buffer = Paso_Coupler_finishCollect(coupler);
1030    
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
     arg.requireWrite();  
1031  #pragma omp parallel for  #pragma omp parallel for
1032      for (dim_t i1 = 0; i1 < m_N1; i1++) {      for (index_t i=0; i<getNumNodes(); i++) {
1033          for (dim_t i0 = 0; i0 < m_N0; i0++) {          const double* src=(m_dofMap[i]<numDOF ?
1034              double* point = arg.getSampleDataRW(i0+m_N0*i1);                  in.getSampleDataRO(m_dofMap[i])
1035              point[0] = xdx.first+i0*xdx.second;                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1036              point[1] = ydy.first+i1*ydy.second;          copy(src, src+numComp, out.getSampleDataRW(i));
         }  
1037      }      }
1038  }  }
1039    
1040  //private  //private
1041  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
1042  {  {
1043      // identifiers are ordered from left to right, bottom to top on each rank,      // identifiers are ordered from left to right, bottom to top globablly.
     // except for the shared nodes which are owned by the rank below / to the  
     // left of the current rank  
1044    
1045      // build node distribution vector first.      // build node distribution vector first.
     // m_nodeDistribution[i] is the first node id on rank i, that is  
1046      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
1047      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1048      m_nodeDistribution[1]=getNumNodes();      const dim_t numDOF=getNumDOF();
1049      for (dim_t k=1; k<m_mpiInfo->size-1; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
1050          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;  
1051      }      }
1052      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
   
1053      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
1054        m_dofId.resize(numDOF);
1055        m_elementId.resize(getNumElements());
1056        m_faceId.resize(getNumFaceElements());
1057    
1058      // the bottom row and left column are not owned by this rank so the  #pragma omp parallel
1059      // identifiers need to be computed accordingly      {
1060      const index_t left = (m_offset0==0 ? 0 : 1);          // nodes
1061      const index_t bottom = (m_offset1==0 ? 0 : 1);  #pragma omp for nowait
1062      if (left>0) {          for (dim_t i1=0; i1<m_N1; i1++) {
1063          const int neighbour=m_mpiInfo->rank-1;              for (dim_t i0=0; i0<m_N0; i0++) {
1064          const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);                  m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;
1065  #pragma omp parallel for              }
         for (dim_t i1=bottom; i1<m_N1; i1++) {  
             m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]  
                 + (i1-bottom+1)*leftN0-1;  
         }  
     }  
     if (bottom>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX;  
         const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);  
 #pragma omp parallel for  
         for (dim_t i0=left; i0<m_N0; i0++) {  
             m_nodeId[i0]=m_nodeDistribution[neighbour]  
                 + (bottomN1-1)*bottomN0 + i0 - left;  
1066          }          }
     }  
     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;  
     }  
1067    
1068      // the rest of the id's are contiguous          // degrees of freedom
1069      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];  #pragma omp for nowait
1070  #pragma omp parallel for          for (dim_t k=0; k<numDOF; k++)
1071      for (dim_t i1=bottom; i1<m_N1; i1++) {              m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;
1072          for (dim_t i0=left; i0<m_N0; i0++) {  
1073              m_nodeId[i0+i1*m_N0] = firstId+i0-left+(i1-bottom)*(m_N0-left);          // elements
1074    #pragma omp for nowait
1075            for (dim_t i1=0; i1<m_NE1; i1++) {
1076                for (dim_t i0=0; i0<m_NE0; i0++) {
1077                    m_elementId[i0+i1*m_NE0]=(m_offset1+i1)*m_gNE0+m_offset0+i0;
1078                }
1079          }          }
     }  
1080    
1081      // element id's          // face elements
1082      m_elementId.resize(getNumElements());  #pragma omp for
1083  #pragma omp parallel for          for (dim_t k=0; k<getNumFaceElements(); k++)
1084      for (dim_t k=0; k<getNumElements(); k++) {              m_faceId[k]=k;
1085          m_elementId[k]=k;      } // end parallel section
     }  
1086    
1087      // face element id's      m_nodeTags.assign(getNumNodes(), 0);
1088      m_faceId.resize(getNumFaceElements());      updateTagsInUse(Nodes);
1089  #pragma omp parallel for  
1090      for (dim_t k=0; k<getNumFaceElements(); k++) {      m_elementTags.assign(getNumElements(), 0);
1091          m_faceId[k]=k;      updateTagsInUse(Elements);
     }  
1092    
1093      // generate face offset vector and set face tags      // generate face offset vector and set face tags
1094      const IndexVector facesPerEdge = getNumFacesPerBoundary();      const IndexVector facesPerEdge = getNumFacesPerBoundary();
# Line 1156  void Rectangle::populateSampleIds() Line 1112  void Rectangle::populateSampleIds()
1112  }  }
1113    
1114  //private  //private
1115  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const  void Rectangle::createPattern()
1116  {  {
1117      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1118      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1119      const int x=node%myN0;      const index_t left = (m_offset0==0 ? 0 : 1);
1120      const int y=node/myN0;      const index_t bottom = (m_offset1==0 ? 0 : 1);
1121      int num=0;  
1122      if (y>0) {      // populate node->DOF mapping with own degrees of freedom.
1123          if (x>0) {      // The rest is assigned in the loop further down
1124              // bottom-left      m_dofMap.assign(getNumNodes(), 0);
1125              index.push_back(node-myN0-1);  #pragma omp parallel for
1126              num++;      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1127          }          for (index_t j=left; j<left+nDOF0; j++) {
1128          // bottom              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
         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++;  
1129          }          }
1130      }      }
1131      if (y<myN1-1) {  
1132          // top      // build list of shared components and neighbours by looping through
1133          index.push_back(node+myN0);      // all potential neighbouring ranks and checking if positions are
1134          num++;      // within bounds
1135          if (x>0) {      const dim_t numDOF=nDOF0*nDOF1;
1136              // top-left      vector<IndexVector> colIndices(numDOF); // for the couple blocks
1137              index.push_back(node+myN0-1);      RankVector neighbour;
1138              num++;      IndexVector offsetInShared(1,0);
1139        IndexVector sendShared, recvShared;
1140        int numShared=0;
1141        const int x=m_mpiInfo->rank%m_NX;
1142        const int y=m_mpiInfo->rank/m_NX;
1143        for (int i1=-1; i1<2; i1++) {
1144            for (int i0=-1; i0<2; i0++) {
1145                // skip this rank
1146                if (i0==0 && i1==0)
1147                    continue;
1148                // location of neighbour rank
1149                const int nx=x+i0;
1150                const int ny=y+i1;
1151                if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {
1152                    neighbour.push_back(ny*m_NX+nx);
1153                    if (i0==0) {
1154                        // sharing top or bottom edge
1155                        const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1156                        const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);
1157                        offsetInShared.push_back(offsetInShared.back()+nDOF0);
1158                        for (dim_t i=0; i<nDOF0; i++, numShared++) {
1159                            sendShared.push_back(firstDOF+i);
1160                            recvShared.push_back(numDOF+numShared);
1161                            if (i>0)
1162                                colIndices[firstDOF+i-1].push_back(numShared);
1163                            colIndices[firstDOF+i].push_back(numShared);
1164                            if (i<nDOF0-1)
1165                                colIndices[firstDOF+i+1].push_back(numShared);
1166                            m_dofMap[firstNode+i]=numDOF+numShared;
1167                        }
1168                    } else if (i1==0) {
1169                        // sharing left or right edge
1170                        const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1171                        const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);
1172                        offsetInShared.push_back(offsetInShared.back()+nDOF1);
1173                        for (dim_t i=0; i<nDOF1; i++, numShared++) {
1174                            sendShared.push_back(firstDOF+i*nDOF0);
1175                            recvShared.push_back(numDOF+numShared);
1176                            if (i>0)
1177                                colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);
1178                            colIndices[firstDOF+i*nDOF0].push_back(numShared);
1179                            if (i<nDOF1-1)
1180                                colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1181                            m_dofMap[firstNode+i*m_N0]=numDOF+numShared;
1182                        }
1183                    } else {
1184                        // sharing a node
1185                        const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1186                        const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);
1187                        offsetInShared.push_back(offsetInShared.back()+1);
1188                        sendShared.push_back(dof);
1189                        recvShared.push_back(numDOF+numShared);
1190                        colIndices[dof].push_back(numShared);
1191                        m_dofMap[node]=numDOF+numShared;
1192                        ++numShared;
1193                    }
1194                }
1195          }          }
1196      }      }
     if (x>0) {  
         // left  
         index.push_back(node-1);  
         num++;  
     }  
1197    
1198      return num;      // create connector
1199  }      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1200                numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
1201                &offsetInShared[0], 1, 0, m_mpiInfo);
1202        Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
1203                numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
1204                &offsetInShared[0], 1, 0, m_mpiInfo);
1205        m_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1206        Paso_SharedComponents_free(snd_shcomp);
1207        Paso_SharedComponents_free(rcv_shcomp);
1208    
1209  //private      // create main and couple blocks
1210  void Rectangle::generateCouplePatterns(Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const      Paso_Pattern *mainPattern = createMainPattern();
1211  {      Paso_Pattern *colPattern, *rowPattern;
1212      IndexVector ptr(1,0);      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
     IndexVector index;  
     const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);  
     const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);  
     const IndexVector faces=getNumFacesPerBoundary();  
1213    
1214      // bottom edge      // allocate paso distribution
1215      dim_t n=0;      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
1216      if (faces[0] == 0) {              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
         index.push_back(2*(myN0+myN1+1));  
         index.push_back(2*(myN0+myN1+1)+1);  
         n+=2;  
         if (faces[2] == 0) {  
             index.push_back(0);  
             index.push_back(1);  
             index.push_back(2);  
             n+=3;  
         }  
     } else if (faces[2] == 0) {  
         index.push_back(1);  
         index.push_back(2);  
         n+=2;  
     }  
     // n=neighbours of bottom-left corner node  
     ptr.push_back(ptr.back()+n);  
     n=0;  
     if (faces[2] == 0) {  
         for (dim_t i=1; i<myN0-1; i++) {  
             index.push_back(i);  
             index.push_back(i+1);  
             index.push_back(i+2);  
             ptr.push_back(ptr.back()+3);  
         }  
         index.push_back(myN0-1);  
         index.push_back(myN0);  
         n+=2;  
         if (faces[1] == 0) {  
             index.push_back(myN0+1);  
             index.push_back(myN0+2);  
             index.push_back(myN0+3);  
             n+=3;  
         }  
     } else {  
         for (dim_t i=1; i<myN0-1; i++) {  
             ptr.push_back(ptr.back());  
         }  
         if (faces[1] == 0) {  
             index.push_back(myN0+2);  
             index.push_back(myN0+3);  
             n+=2;  
         }  
     }  
     // n=neighbours of bottom-right corner node  
     ptr.push_back(ptr.back()+n);  
1217    
1218      // 2nd row to 2nd last row      // finally create the system matrix
1219      for (dim_t i=1; i<myN1-1; i++) {      m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
1220          // left edge              distribution, distribution, mainPattern, colPattern, rowPattern,
1221          if (faces[0] == 0) {              m_connector, m_connector);
             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());  
         }  
     }  
1222    
1223      // top edge      Paso_Distribution_free(distribution);
     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);  
1224    
1225      dim_t M=ptr.size()-1;      // useful debug output
1226      map<index_t,index_t> idMap;      /*
1227      dim_t N=0;      cout << "--- rcv_shcomp ---" << endl;
1228      for (IndexVector::iterator it=index.begin(); it!=index.end(); it++) {      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
1229          if (idMap.find(*it)==idMap.end()) {      for (size_t i=0; i<neighbour.size(); i++) {
1230              idMap[*it]=N;          cout << "neighbor[" << i << "]=" << neighbour[i]
1231              *it=N++;              << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
         } else {  
             *it=idMap[*it];  
         }  
1232      }      }
1233        for (size_t i=0; i<recvShared.size(); i++) {
1234            cout << "shared[" << i << "]=" << recvShared[i] << endl;
1235        }
1236        cout << "--- snd_shcomp ---" << endl;
1237        for (size_t i=0; i<sendShared.size(); i++) {
1238            cout << "shared[" << i << "]=" << sendShared[i] << endl;
1239        }
1240        cout << "--- dofMap ---" << endl;
1241        for (size_t i=0; i<m_dofMap.size(); i++) {
1242            cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
1243        }
1244        cout << "--- colIndices ---" << endl;
1245        for (size_t i=0; i<colIndices.size(); i++) {
1246            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
1247        }
1248        */
1249    
1250      /*      /*
1251      cout << "--- colCouple_pattern ---" << endl;      cout << "--- main_pattern ---" << endl;
1252      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
1253      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<mainPattern->numOutput+1; i++) {
1254          cout << "ptr[" << i << "]=" << ptr[i] << endl;          cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
1255      }      }
1256      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
1257          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
1258      }      }
1259      */      */
1260    
1261      // now build the row couple pattern      /*
1262      IndexVector ptr2(1,0);      cout << "--- colCouple_pattern ---" << endl;
1263      IndexVector index2;      cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
1264      for (dim_t id=0; id<N; id++) {      for (size_t i=0; i<colPattern->numOutput+1; i++) {
1265          n=0;          cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
1266          for (dim_t i=0; i<M; i++) {      }
1267              for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {      for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
1268                  if (index[j]==id) {          cout << "index[" << i << "]=" << colPattern->index[i] << endl;
                     index2.push_back(i);  
                     n++;  
                     break;  
                 }  
             }  
         }  
         ptr2.push_back(ptr2.back()+n);  
1269      }      }
1270        */
1271    
1272      /*      /*
1273      cout << "--- rowCouple_pattern ---" << endl;      cout << "--- rowCouple_pattern ---" << endl;
1274      cout << "M=" << N << ", N=" << M << endl;      cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
1275      for (size_t i=0; i<ptr2.size(); i++) {      for (size_t i=0; i<rowPattern->numOutput+1; i++) {
1276          cout << "ptr[" << i << "]=" << ptr2[i] << endl;          cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
1277      }      }
1278      for (size_t i=0; i<index2.size(); i++) {      for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
1279          cout << "index[" << i << "]=" << index2[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
1280      }      }
1281      */      */
1282    
1283      // paso will manage the memory      Paso_Pattern_free(mainPattern);
1284      index_t* indexC = MEMALLOC(index.size(), index_t);      Paso_Pattern_free(colPattern);
1285      index_t* ptrC = MEMALLOC(ptr.size(), index_t);      Paso_Pattern_free(rowPattern);
     copy(index.begin(), index.end(), indexC);  
     copy(ptr.begin(), ptr.end(), ptrC);  
     *colPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);  
   
     // paso will manage the memory  
     indexC = MEMALLOC(index2.size(), index_t);  
     ptrC = MEMALLOC(ptr2.size(), index_t);  
     copy(index2.begin(), index2.end(), indexC);  
     copy(ptr2.begin(), ptr2.end(), ptrC);  
     *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);  
1286  }  }
1287    
1288  //protected  //protected
# Line 1411  void Rectangle::interpolateNodesOnElemen Line 1291  void Rectangle::interpolateNodesOnElemen
1291  {  {
1292      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1293      if (reduced) {      if (reduced) {
1294          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */          out.requireWrite();
1295          const double tmp0_0 = 0.25000000000000000000;          const double c0 = .25;
1296  #pragma omp parallel for  #pragma omp parallel for
1297          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1298              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 1422  void Rectangle::interpolateNodesOnElemen Line 1302  void Rectangle::interpolateNodesOnElemen
1302                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1303                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1304                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1305                      o[INDEX2(i,numComp,0)] = tmp0_0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1306                  } /* end of component loop i */                  } /* end of component loop i */
1307              } /* end of k0 loop */              } /* end of k0 loop */
1308          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */  
1309      } else {      } else {
1310          /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */          out.requireWrite();
1311          const double tmp0_2 = 0.62200846792814621559;          const double c0 = .16666666666666666667;
1312          const double tmp0_1 = 0.044658198738520451079;          const double c1 = .044658198738520451079;
1313          const double tmp0_0 = 0.16666666666666666667;          const double c2 = .62200846792814621559;
1314  #pragma omp parallel for  #pragma omp parallel for
1315          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1316              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 1441  void Rectangle::interpolateNodesOnElemen Line 1320  void Rectangle::interpolateNodesOnElemen
1320                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1321                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1322                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1323                      o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);
1324                      o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);                      o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);
1325                      o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);                      o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);
1326                      o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);                      o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);
1327                  } /* end of component loop i */                  } /* end of component loop i */
1328              } /* end of k0 loop */              } /* end of k0 loop */
1329          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */  
1330      }      }
1331  }  }
1332    
# Line 1458  void Rectangle::interpolateNodesOnFaces( Line 1336  void Rectangle::interpolateNodesOnFaces(
1336  {  {
1337      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1338      if (reduced) {      if (reduced) {
1339          /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */          out.requireWrite();
1340          if (m_faceOffset[0] > -1) {          const double c0 = .5;
1341              const double tmp0_0 = 0.50000000000000000000;  #pragma omp parallel
1342  #pragma omp parallel for          {
1343              for (index_t k1=0; k1 < m_NE1; ++k1) {              if (m_faceOffset[0] > -1) {
1344                  const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));  #pragma omp for nowait
1345                  const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1346                  double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1347                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1348                      o[INDEX2(i,numComp,0)] = tmp0_0*(f_00[i] + f_01[i]);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1349                  } /* end of component loop i */                      for (index_t i=0; i < numComp; ++i) {
1350              } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
1351          } /* end of face 0 */                      } /* end of component loop i */
1352          if (m_faceOffset[1] > -1) {                  } /* end of k1 loop */
1353              const double tmp0_0 = 0.50000000000000000000;              } /* end of face 0 */
1354  #pragma omp parallel for              if (m_faceOffset[1] > -1) {
1355              for (index_t k1=0; k1 < m_NE1; ++k1) {  #pragma omp for nowait
1356                  const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1357                  const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1358                  double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1359                  for (index_t i=0; i < numComp; ++i) {                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1360                      o[INDEX2(i,numComp,0)] = tmp0_0*(f_10[i] + f_11[i]);                      for (index_t i=0; i < numComp; ++i) {
1361                  } /* end of component loop i */                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
1362              } /* end of k1 loop */                      } /* end of component loop i */
1363          } /* end of face 1 */                  } /* end of k1 loop */
1364          if (m_faceOffset[2] > -1) {              } /* end of face 1 */
1365              const double tmp0_0 = 0.50000000000000000000;              if (m_faceOffset[2] > -1) {
1366  #pragma omp parallel for  #pragma omp for nowait
1367              for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1368                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1369                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1370                  double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1371                  for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1372                      o[INDEX2(i,numComp,0)] = tmp0_0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
1373                  } /* end of component loop i */                      } /* end of component loop i */
1374              } /* end of k0 loop */                  } /* end of k0 loop */
1375          } /* end of face 2 */              } /* end of face 2 */
1376          if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1377              const double tmp0_0 = 0.50000000000000000000;  #pragma omp for nowait
1378  #pragma omp parallel for                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1379              for (index_t k0=0; k0 < m_NE0; ++k0) {                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1380                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1381                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1382                  double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      for (index_t i=0; i < numComp; ++i) {
1383                  for (index_t i=0; i < numComp; ++i) {                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1384                      o[INDEX2(i,numComp,0)] = tmp0_0*(f_01[i] + f_11[i]);                      } /* end of component loop i */
1385                  } /* end of component loop i */                  } /* end of k0 loop */
1386              } /* end of k0 loop */              } /* end of face 3 */
1387          } /* end of face 3 */          } // end of parallel section
         /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */  
1388      } else {      } else {
1389          /* GENERATOR SNIP_INTERPOLATE_FACES TOP */          out.requireWrite();
1390          if (m_faceOffset[0] > -1) {          const double c0 = 0.21132486540518711775;
1391              const double tmp0_1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1392              const double tmp0_0 = 0.21132486540518711775;  #pragma omp parallel
1393  #pragma omp parallel for          {
1394              for (index_t k1=0; k1 < m_NE1; ++k1) {              if (m_faceOffset[0] > -1) {
1395                  const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));  #pragma omp for nowait
1396                  const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1397                  double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1398                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1399                      o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_01[i]*tmp0_0;                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1400                      o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_01[i]*tmp0_1;                      for (index_t i=0; i < numComp; ++i) {
1401                  } /* end of component loop i */                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;
1402              } /* end of k1 loop */                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_01[i]*c1;
1403          } /* end of face 0 */                      } /* end of component loop i */
1404          if (m_faceOffset[1] > -1) {                  } /* end of k1 loop */
1405              const double tmp0_1 = 0.21132486540518711775;              } /* end of face 0 */
1406              const double tmp0_0 = 0.78867513459481288225;              if (m_faceOffset[1] > -1) {
1407  #pragma omp parallel for  #pragma omp for nowait
1408              for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1409                  const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1410                  const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1411                  double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1412                  for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1413                      o[INDEX2(i,numComp,0)] = f_10[i]*tmp0_0 + f_11[i]*tmp0_1;                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;
1414                      o[INDEX2(i,numComp,1)] = f_10[i]*tmp0_1 + f_11[i]*tmp0_0;                          o[INDEX2(i,numComp,1)] = f_10[i]*c0 + f_11[i]*c1;
1415                  } /* end of component loop i */                      } /* end of component loop i */
1416              } /* end of k1 loop */                  } /* end of k1 loop */
1417          } /* end of face 1 */              } /* end of face 1 */
1418          if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1419              const double tmp0_1 = 0.78867513459481288225;  #pragma omp for nowait
1420              const double tmp0_0 = 0.21132486540518711775;                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1421  #pragma omp parallel for                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1422              for (index_t k0=0; k0 < m_NE0; ++k0) {                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1423                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1424                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      for (index_t i=0; i < numComp; ++i) {
1425                  double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;
1426                  for (index_t i=0; i < numComp; ++i) {                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_10[i]*c1;
1427                      o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_10[i]*tmp0_0;                      } /* end of component loop i */
1428                      o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;                  } /* end of k0 loop */
1429                  } /* end of component loop i */              } /* end of face 2 */
1430              } /* end of k0 loop */              if (m_faceOffset[3] > -1) {
1431          } /* end of face 2 */  #pragma omp for nowait
1432          if (m_faceOffset[3] > -1) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1433              const double tmp0_1 = 0.78867513459481288225;                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1434              const double tmp0_0 = 0.21132486540518711775;                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1435  #pragma omp parallel for                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1436              for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t i=0; i < numComp; ++i) {
1437                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;
1438                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                          o[INDEX2(i,numComp,1)] = f_01[i]*c0 + f_11[i]*c1;
1439                  double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      } /* end of component loop i */
1440                  for (index_t i=0; i < numComp; ++i) {                  } /* end of k0 loop */
1441                      o[INDEX2(i,numComp,0)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;              } /* end of face 3 */
1442                      o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_0 + f_11[i]*tmp0_1;          } // end of parallel section
                 } /* end of component loop i */  
             } /* end of k0 loop */  
         } /* end of face 3 */  
         /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */  
1443      }      }
1444  }  }
1445    
1446    //protected
1447    void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1448            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1449            const escript::Data& C, const escript::Data& D,
1450            const escript::Data& X, const escript::Data& Y,
1451            const escript::Data& d, const escript::Data& y) const
1452    {
1453        const double h0 = m_l0/m_gNE0;
1454        const double h1 = m_l1/m_gNE1;
1455        const double w0 = -0.1555021169820365539*h1/h0;
1456        const double w1 = 0.041666666666666666667;
1457        const double w2 = -0.15550211698203655390;
1458        const double w3 = 0.041666666666666666667*h0/h1;
1459        const double w4 = 0.15550211698203655390;
1460        const double w5 = -0.041666666666666666667;
1461        const double w6 = -0.01116454968463011277*h1/h0;
1462        const double w7 = 0.011164549684630112770;
1463        const double w8 = -0.011164549684630112770;
1464        const double w9 = -0.041666666666666666667*h1/h0;
1465        const double w10 = -0.041666666666666666667*h0/h1;
1466        const double w11 = 0.1555021169820365539*h1/h0;
1467        const double w12 = 0.1555021169820365539*h0/h1;
1468        const double w13 = 0.01116454968463011277*h0/h1;
1469        const double w14 = 0.01116454968463011277*h1/h0;
1470        const double w15 = 0.041666666666666666667*h1/h0;
1471        const double w16 = -0.01116454968463011277*h0/h1;
1472        const double w17 = -0.1555021169820365539*h0/h1;
1473        const double w18 = -0.33333333333333333333*h1/h0;
1474        const double w19 = 0.25000000000000000000;
1475        const double w20 = -0.25000000000000000000;
1476        const double w21 = 0.16666666666666666667*h0/h1;
1477        const double w22 = -0.16666666666666666667*h1/h0;
1478        const double w23 = -0.16666666666666666667*h0/h1;
1479        const double w24 = 0.33333333333333333333*h1/h0;
1480        const double w25 = 0.33333333333333333333*h0/h1;
1481        const double w26 = 0.16666666666666666667*h1/h0;
1482        const double w27 = -0.33333333333333333333*h0/h1;
1483        const double w28 = -0.032861463941450536761*h1;
1484        const double w29 = -0.032861463941450536761*h0;
1485        const double w30 = -0.12264065304058601714*h1;
1486        const double w31 = -0.0023593469594139828636*h1;
1487        const double w32 = -0.008805202725216129906*h0;
1488        const double w33 = -0.008805202725216129906*h1;
1489        const double w34 = 0.032861463941450536761*h1;
1490        const double w35 = 0.008805202725216129906*h1;
1491        const double w36 = 0.008805202725216129906*h0;
1492        const double w37 = 0.0023593469594139828636*h1;
1493        const double w38 = 0.12264065304058601714*h1;
1494        const double w39 = 0.032861463941450536761*h0;
1495        const double w40 = -0.12264065304058601714*h0;
1496        const double w41 = -0.0023593469594139828636*h0;
1497        const double w42 = 0.0023593469594139828636*h0;
1498        const double w43 = 0.12264065304058601714*h0;
1499        const double w44 = -0.16666666666666666667*h1;
1500        const double w45 = -0.083333333333333333333*h0;
1501        const double w46 = 0.083333333333333333333*h1;
1502        const double w47 = 0.16666666666666666667*h1;
1503        const double w48 = 0.083333333333333333333*h0;
1504        const double w49 = -0.16666666666666666667*h0;
1505        const double w50 = 0.16666666666666666667*h0;
1506        const double w51 = -0.083333333333333333333*h1;
1507        const double w52 = 0.025917019497006092316*h0*h1;
1508        const double w53 = 0.0018607582807716854616*h0*h1;
1509        const double w54 = 0.0069444444444444444444*h0*h1;
1510        const double w55 = 0.09672363354357992482*h0*h1;
1511        const double w56 = 0.00049858867864229740201*h0*h1;
1512        const double w57 = 0.055555555555555555556*h0*h1;
1513        const double w58 = 0.027777777777777777778*h0*h1;
1514        const double w59 = 0.11111111111111111111*h0*h1;
1515        const double w60 = -0.19716878364870322056*h1;
1516        const double w61 = -0.19716878364870322056*h0;
1517        const double w62 = -0.052831216351296779436*h0;
1518        const double w63 = -0.052831216351296779436*h1;
1519        const double w64 = 0.19716878364870322056*h1;
1520        const double w65 = 0.052831216351296779436*h1;
1521        const double w66 = 0.19716878364870322056*h0;
1522        const double w67 = 0.052831216351296779436*h0;
1523        const double w68 = -0.5*h1;
1524        const double w69 = -0.5*h0;
1525        const double w70 = 0.5*h1;
1526        const double w71 = 0.5*h0;
1527        const double w72 = 0.1555021169820365539*h0*h1;
1528        const double w73 = 0.041666666666666666667*h0*h1;
1529        const double w74 = 0.01116454968463011277*h0*h1;
1530        const double w75 = 0.25*h0*h1;
1531    
1532        rhs.requireWrite();
1533    #pragma omp parallel
1534        {
1535            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1536    #pragma omp for
1537                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
1538                    for (index_t k0=0; k0<m_NE0; ++k0)  {
1539                        bool add_EM_S=false;
1540                        bool add_EM_F=false;
1541                        vector<double> EM_S(4*4, 0);
1542                        vector<double> EM_F(4, 0);
1543                        const index_t e = k0 + m_NE0*k1;
1544                        ///////////////
1545                        // process A //
1546                        ///////////////
1547                        if (!A.isEmpty()) {
1548                            add_EM_S=true;
1549                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1550                            if (A.actsExpanded()) {
1551                                const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1552                                const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1553                                const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1554                                const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1555                                const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1556                                const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1557                                const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1558                                const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1559                                const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1560                                const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1561                                const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1562                                const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1563                                const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1564                                const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1565                                const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1566                                const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1567                                const register double tmp0_0 = A_01_0 + A_01_3;
1568                                const register double tmp1_0 = A_00_0 + A_00_1;
1569                                const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1570                                const register double tmp3_0 = A_00_2 + A_00_3;
1571                                const register double tmp4_0 = A_10_1 + A_10_2;
1572                                const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1573                                const register double tmp6_0 = A_01_3 + A_10_0;
1574                                const register double tmp7_0 = A_01_0 + A_10_3;
1575                                const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1576                                const register double tmp9_0 = A_01_0 + A_10_0;
1577                                const register double tmp12_0 = A_11_0 + A_11_2;
1578                                const register double tmp10_0 = A_01_3 + A_10_3;
1579                                const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
1580                                const register double tmp13_0 = A_01_2 + A_10_1;
1581                                const register double tmp11_0 = A_11_1 + A_11_3;
1582                                const register double tmp18_0 = A_01_1 + A_10_1;
1583                                const register double tmp15_0 = A_01_1 + A_10_2;
1584                                const register double tmp16_0 = A_10_0 + A_10_3;
1585                                const register double tmp17_0 = A_01_1 + A_01_2;
1586                                const register double tmp19_0 = A_01_2 + A_10_2;
1587                                const register double tmp0_1 = A_10_3*w8;
1588                                const register double tmp1_1 = tmp0_0*w1;
1589                                const register double tmp2_1 = A_01_1*w4;
1590                                const register double tmp3_1 = tmp1_0*w0;
1591                                const register double tmp4_1 = A_01_2*w7;
1592                                const register double tmp5_1 = tmp2_0*w3;
1593                                const register double tmp6_1 = tmp3_0*w6;
1594                                const register double tmp7_1 = A_10_0*w2;
1595                                const register double tmp8_1 = tmp4_0*w5;
1596                                const register double tmp9_1 = tmp2_0*w10;
1597                                const register double tmp14_1 = A_10_0*w8;
1598                                const register double tmp23_1 = tmp3_0*w14;
1599                                const register double tmp35_1 = A_01_0*w8;
1600                                const register double tmp54_1 = tmp13_0*w8;
1601                                const register double tmp20_1 = tmp9_0*w4;
1602                                const register double tmp25_1 = tmp12_0*w12;
1603                                const register double tmp44_1 = tmp7_0*w7;
1604                                const register double tmp26_1 = tmp10_0*w4;
1605                                const register double tmp52_1 = tmp18_0*w8;
1606                                const register double tmp48_1 = A_10_1*w7;
1607                                const register double tmp46_1 = A_01_3*w8;
1608                                const register double tmp50_1 = A_01_0*w2;
1609                                const register double tmp56_1 = tmp19_0*w8;
1610                                const register double tmp19_1 = A_10_3*w2;
1611                                const register double tmp47_1 = A_10_2*w4;
1612                                const register double tmp16_1 = tmp3_0*w0;
1613                                const register double tmp18_1 = tmp1_0*w6;
1614                                const register double tmp31_1 = tmp11_0*w12;
1615                                const register double tmp55_1 = tmp15_0*w2;
1616                                const register double tmp39_1 = A_10_2*w7;
1617                                const register double tmp11_1 = tmp6_0*w7;
1618                                const register double tmp40_1 = tmp11_0*w17;
1619                                const register double tmp34_1 = tmp15_0*w8;
1620                                const register double tmp33_1 = tmp14_0*w5;
1621                                const register double tmp24_1 = tmp11_0*w13;
1622                                const register double tmp43_1 = tmp17_0*w5;
1623                                const register double tmp15_1 = A_01_2*w4;
1624                                const register double tmp53_1 = tmp19_0*w2;
1625                                const register double tmp27_1 = tmp3_0*w11;
1626                                const register double tmp32_1 = tmp13_0*w2;
1627                                const register double tmp10_1 = tmp5_0*w9;
1628                                const register double tmp37_1 = A_10_1*w4;
1629                                const register double tmp38_1 = tmp5_0*w15;
1630                                const register double tmp17_1 = A_01_1*w7;
1631                                const register double tmp12_1 = tmp7_0*w4;
1632                                const register double tmp22_1 = tmp10_0*w7;
1633                                const register double tmp57_1 = tmp18_0*w2;
1634                                const register double tmp28_1 = tmp9_0*w7;
1635                                const register double tmp29_1 = tmp1_0*w14;
1636                                const register double tmp51_1 = tmp11_0*w16;
1637                                const register double tmp42_1 = tmp12_0*w16;
1638                                const register double tmp49_1 = tmp12_0*w17;
1639                                const register double tmp21_1 = tmp1_0*w11;
1640                                const register double tmp45_1 = tmp6_0*w4;
1641                                const register double tmp13_1 = tmp8_0*w1;
1642                                const register double tmp36_1 = tmp16_0*w1;
1643                                const register double tmp41_1 = A_01_3*w2;
1644                                const register double tmp30_1 = tmp12_0*w13;
1645                                EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1646                                EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
1647                                EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
1648                                EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
1649                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
1650                                EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
1651                                EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
1652                                EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
1653                                EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1654                                EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
1655                                EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
1656                                EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
1657                                EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
1658                                EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
1659                                EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
1660                                EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1661                            } else { // constant data
1662                                const register double A_00 = A_p[INDEX2(0,0,2)];
1663                                const register double A_10 = A_p[INDEX2(1,0,2)];
1664                                const register double A_01 = A_p[INDEX2(0,1,2)];
1665                                const register double A_11 = A_p[INDEX2(1,1,2)];
1666                                const register double tmp0_0 = A_01 + A_10;
1667                                const register double tmp0_1 = A_00*w18;
1668                                const register double tmp1_1 = A_01*w19;
1669                                const register double tmp2_1 = A_10*w20;
1670                                const register double tmp3_1 = A_11*w21;
1671                                const register double tmp4_1 = A_00*w22;
1672                                const register double tmp5_1 = tmp0_0*w19;
1673                                const register double tmp6_1 = A_11*w23;
1674                                const register double tmp7_1 = A_11*w25;
1675                                const register double tmp8_1 = A_00*w24;
1676                                const register double tmp9_1 = tmp0_0*w20;
1677                                const register double tmp10_1 = A_01*w20;
1678                                const register double tmp11_1 = A_11*w27;
1679                                const register double tmp12_1 = A_00*w26;
1680                                const register double tmp13_1 = A_10*w19;
1681                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1682                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1683                                EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1684                                EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1685                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1686                                EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1687                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1688                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1689                                EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1690                                EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1691                                EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1692                                EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1693                                EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1694                                EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1695                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1696                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1697                            }
1698                        }
1699                        ///////////////
1700                        // process B //
1701                        ///////////////
1702                        if (!B.isEmpty()) {
1703                            add_EM_S=true;
1704                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1705                            if (B.actsExpanded()) {
1706                                const register double B_0_0 = B_p[INDEX2(0,0,2)];
1707                                const register double B_1_0 = B_p[INDEX2(1,0,2)];
1708                                const register double B_0_1 = B_p[INDEX2(0,1,2)];
1709                                const register double B_1_1 = B_p[INDEX2(1,1,2)];
1710                                const register double B_0_2 = B_p[INDEX2(0,2,2)];
1711                                const register double B_1_2 = B_p[INDEX2(1,2,2)];
1712                                const register double B_0_3 = B_p[INDEX2(0,3,2)];
1713                                const register double B_1_3 = B_p[INDEX2(1,3,2)];
1714                                const register double tmp0_0 = B_1_0 + B_1_1;
1715                                const register double tmp1_0 = B_1_2 + B_1_3;
1716                                const register double tmp2_0 = B_0_1 + B_0_3;
1717                                const register double tmp3_0 = B_0_0 + B_0_2;
1718                                const register double tmp63_1 = B_1_1*w42;
1719                                const register double tmp79_1 = B_1_1*w40;
1720                                const register double tmp37_1 = tmp3_0*w35;
1721                                const register double tmp8_1 = tmp0_0*w32;
1722                                const register double tmp71_1 = B_0_1*w34;
1723                                const register double tmp19_1 = B_0_3*w31;
1724                                const register double tmp15_1 = B_0_3*w34;
1725                                const register double tmp9_1 = tmp3_0*w34;
1726                                const register double tmp35_1 = B_1_0*w36;
1727                                const register double tmp66_1 = B_0_3*w28;
1728                                const register double tmp28_1 = B_1_0*w42;
1729                                const register double tmp22_1 = B_1_0*w40;
1730                                const register double tmp16_1 = B_1_2*w29;
1731                                const register double tmp6_1 = tmp2_0*w35;
1732                                const register double tmp55_1 = B_1_3*w40;
1733                                const register double tmp50_1 = B_1_3*w42;
1734                                const register double tmp7_1 = tmp1_0*w29;
1735                                const register double tmp1_1 = tmp1_0*w32;
1736                                const register double tmp57_1 = B_0_3*w30;
1737                                const register double tmp18_1 = B_1_1*w32;
1738                                const register double tmp53_1 = B_1_0*w41;
1739                                const register double tmp61_1 = B_1_3*w36;
1740                                const register double tmp27_1 = B_0_3*w38;
1741                                const register double tmp64_1 = B_0_2*w30;
1742                                const register double tmp76_1 = B_0_1*w38;
1743                                const register double tmp39_1 = tmp2_0*w34;
1744                                const register double tmp62_1 = B_0_1*w31;
1745                                const register double tmp56_1 = B_0_0*w31;
1746                                const register double tmp49_1 = B_1_1*w36;
1747                                const register double tmp2_1 = B_0_2*w31;
1748                                const register double tmp23_1 = B_0_2*w33;
1749                                const register double tmp38_1 = B_1_1*w43;
1750                                const register double tmp74_1 = B_1_2*w41;
1751                                const register double tmp43_1 = B_1_1*w41;
1752                                const register double tmp58_1 = B_0_2*w28;
1753                                const register double tmp67_1 = B_0_0*w33;
1754                                const register double tmp33_1 = tmp0_0*w39;
1755                                const register double tmp4_1 = B_0_0*w28;
1756                                const register double tmp20_1 = B_0_0*w30;
1757                                const register double tmp13_1 = B_0_2*w38;
1758                                const register double tmp65_1 = B_1_2*w43;
1759                                const register double tmp0_1 = tmp0_0*w29;
1760                                const register double tmp41_1 = tmp3_0*w33;
1761                                const register double tmp73_1 = B_0_2*w37;
1762                                const register double tmp69_1 = B_0_0*w38;
1763                                const register double tmp48_1 = B_1_2*w39;
1764                                const register double tmp59_1 = B_0_1*w33;
1765                                const register double tmp17_1 = B_1_3*w41;
1766                                const register double tmp5_1 = B_0_3*w33;
1767                                const register double tmp3_1 = B_0_1*w30;
1768                                const register double tmp21_1 = B_0_1*w28;
1769                                const register double tmp42_1 = B_1_0*w29;
1770                                const register double tmp54_1 = B_1_2*w32;
1771                                const register double tmp60_1 = B_1_0*w39;
1772                                const register double tmp32_1 = tmp1_0*w36;
1773                                const register double tmp10_1 = B_0_1*w37;
1774                                const register double tmp14_1 = B_0_0*w35;
1775                                const register double tmp29_1 = B_0_1*w35;
1776                                const register double tmp26_1 = B_1_2*w36;
1777                                const register double tmp30_1 = B_1_3*w43;
1778                                const register double tmp70_1 = B_0_2*w35;
1779                                const register double tmp34_1 = B_1_3*w39;
1780                                const register double tmp51_1 = B_1_0*w43;
1781                                const register double tmp31_1 = B_0_2*w34;
1782                                const register double tmp45_1 = tmp3_0*w28;
1783                                const register double tmp11_1 = tmp1_0*w39;
1784                                const register double tmp52_1 = B_1_1*w29;
1785                                const register double tmp44_1 = B_1_3*w32;
1786                                const register double tmp25_1 = B_1_1*w39;
1787                                const register double tmp47_1 = tmp2_0*w33;
1788                                const register double tmp72_1 = B_1_3*w29;
1789                                const register double tmp40_1 = tmp2_0*w28;
1790                                const register double tmp46_1 = B_1_2*w40;
1791                                const register double tmp36_1 = B_1_2*w42;
1792                                const register double tmp24_1 = B_0_0*w37;
1793                                const register double tmp77_1 = B_0_3*w35;
1794                                const register double tmp68_1 = B_0_3*w37;
1795                                const register double tmp78_1 = B_0_0*w34;
1796                                const register double tmp12_1 = tmp0_0*w36;
1797                                const register double tmp75_1 = B_1_0*w32;
1798                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1799                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1800                                EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1801                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
1802                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1803                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1804                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1805                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1806                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1807                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1808                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1809                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1810                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
1811                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1812                                EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1813                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1814                            } else { // constant data
1815                                const register double B_0 = B_p[0];
1816                                const register double B_1 = B_p[1];
1817                                const register double tmp0_1 = B_0*w44;
1818                                const register double tmp1_1 = B_1*w45;
1819                                const register double tmp2_1 = B_0*w46;
1820                                const register double tmp3_1 = B_0*w47;
1821                                const register double tmp4_1 = B_1*w48;
1822                                const register double tmp5_1 = B_1*w49;
1823                                const register double tmp6_1 = B_1*w50;
1824                                const register double tmp7_1 = B_0*w51;
1825                                EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1826                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1827                                EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
1828                                EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
1829                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1830                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
1831                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
1832                                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
1833                                EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
1834                                EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
1835                                EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
1836                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1837                                EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
1838                                EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
1839                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
1840                                EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
1841                            }
1842                        }
1843                        ///////////////
1844                        // process C //
1845                        ///////////////
1846                        if (!C.isEmpty()) {
1847                            add_EM_S=true;
1848                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
1849                            if (C.actsExpanded()) {
1850                                const register double C_0_0 = C_p[INDEX2(0,0,2)];
1851                                const register double C_1_0 = C_p[INDEX2(1,0,2)];
1852                                const register double C_0_1 = C_p[INDEX2(0,1,2)];
1853                                const register double C_1_1 = C_p[INDEX2(1,1,2)];
1854                                const register double C_0_2 = C_p[INDEX2(0,2,2)];
1855                                const register double C_1_2 = C_p[INDEX2(1,2,2)];
1856                                const register double C_0_3 = C_p[INDEX2(0,3,2)];
1857                                const register double C_1_3 = C_p[INDEX2(1,3,2)];
1858                                const register double tmp0_0 = C_1_0 + C_1_1;
1859                                const register double tmp1_0 = C_1_2 + C_1_3;
1860                                const register double tmp2_0 = C_0_1 + C_0_3;
1861                                const register double tmp3_0 = C_0_0 + C_0_2;
1862                                const register double tmp64_1 = C_0_2*w30;
1863                                const register double tmp14_1 = C_0_2*w28;
1864                                const register double tmp19_1 = C_0_3*w31;
1865                                const register double tmp22_1 = C_1_0*w40;
1866                                const register double tmp37_1 = tmp3_0*w35;
1867                                const register double tmp29_1 = C_0_1*w35;
1868                                const register double tmp73_1 = C_0_2*w37;
1869                                const register double tmp74_1 = C_1_2*w41;
1870                                const register double tmp52_1 = C_1_3*w39;
1871                                const register double tmp25_1 = C_1_1*w39;
1872                                const register double tmp62_1 = C_0_1*w31;
1873                                const register double tmp79_1 = C_1_1*w40;
1874                                const register double tmp43_1 = C_1_1*w36;
1875                                const register double tmp27_1 = C_0_3*w38;
1876                                const register double tmp28_1 = C_1_0*w42;
1877                                const register double tmp63_1 = C_1_1*w42;
1878                                const register double tmp59_1 = C_0_3*w34;
1879                                const register double tmp72_1 = C_1_3*w29;
1880                                const register double tmp40_1 = tmp2_0*w35;
1881                                const register double tmp13_1 = C_0_3*w30;
1882                                const register double tmp51_1 = C_1_2*w40;
1883                                const register double tmp54_1 = C_1_2*w42;
1884                                const register double tmp12_1 = C_0_0*w31;
1885                                const register double tmp2_1 = tmp1_0*w32;
1886                                const register double tmp68_1 = C_0_2*w31;
1887                                const register double tmp75_1 = C_1_0*w32;
1888                                const register double tmp49_1 = C_1_1*w41;
1889                                const register double tmp4_1 = C_0_2*w35;
1890                                const register double tmp66_1 = C_0_3*w28;
1891                                const register double tmp56_1 = C_0_1*w37;
1892                                const register double tmp5_1 = C_0_1*w34;
1893                                const register double tmp38_1 = tmp2_0*w34;
1894                                const register double tmp76_1 = C_0_1*w38;
1895                                const register double tmp21_1 = C_0_1*w28;
1896                                const register double tmp69_1 = C_0_1*w30;
1897                                const register double tmp53_1 = C_1_0*w36;
1898                                const register double tmp42_1 = C_1_2*w39;
1899                                const register double tmp32_1 = tmp1_0*w29;
1900                                const register double tmp45_1 = C_1_0*w43;
1901                                const register double tmp33_1 = tmp0_0*w32;
1902                                const register double tmp35_1 = C_1_0*w41;
1903                                const register double tmp26_1 = C_1_2*w36;
1904                                const register double tmp67_1 = C_0_0*w33;
1905                                const register double tmp31_1 = C_0_2*w34;
1906                                const register double tmp20_1 = C_0_0*w30;
1907                                const register double tmp70_1 = C_0_0*w28;
1908                                const register double tmp8_1 = tmp0_0*w39;
1909                                const register double tmp30_1 = C_1_3*w43;
1910                                const register double tmp0_1 = tmp0_0*w29;
1911                                const register double tmp17_1 = C_1_3*w41;
1912                                const register double tmp58_1 = C_0_0*w35;
1913                                const register double tmp9_1 = tmp3_0*w33;
1914                                const register double tmp61_1 = C_1_3*w36;
1915                                const register double tmp41_1 = tmp3_0*w34;
1916                                const register double tmp50_1 = C_1_3*w32;
1917                                const register double tmp18_1 = C_1_1*w32;
1918                                const register double tmp6_1 = tmp1_0*w36;
1919                                const register double tmp3_1 = C_0_0*w38;
1920                                const register double tmp34_1 = C_1_1*w29;
1921                                const register double tmp77_1 = C_0_3*w35;
1922                                const register double tmp65_1 = C_1_2*w43;
1923                                const register double tmp71_1 = C_0_3*w33;
1924                                const register double tmp55_1 = C_1_1*w43;
1925                                const register double tmp46_1 = tmp3_0*w28;
1926                                const register double tmp24_1 = C_0_0*w37;
1927                                const register double tmp10_1 = tmp1_0*w39;
1928                                const register double tmp48_1 = C_1_0*w29;
1929                                const register double tmp15_1 = C_0_1*w33;
1930                                const register double tmp36_1 = C_1_2*w32;
1931                                const register double tmp60_1 = C_1_0*w39;
1932                                const register double tmp47_1 = tmp2_0*w33;
1933                                const register double tmp16_1 = C_1_2*w29;
1934                                const register double tmp1_1 = C_0_3*w37;
1935                                const register double tmp7_1 = tmp2_0*w28;
1936                                const register double tmp39_1 = C_1_3*w40;
1937                                const register double tmp44_1 = C_1_3*w42;
1938                                const register double tmp57_1 = C_0_2*w38;
1939                                const register double tmp78_1 = C_0_0*w34;
1940                                const register double tmp11_1 = tmp0_0*w36;
1941                                const register double tmp23_1 = C_0_2*w33;
1942                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1943                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1944                                EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1945                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
1946                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1947                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1948                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1949                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1950                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1951                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1952                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1953                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1954                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
1955                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1956                                EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1957                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1958                            } else { // constant data
1959                                const register double C_0 = C_p[0];
1960                                const register double C_1 = C_p[1];
1961                                const register double tmp0_1 = C_0*w47;
1962                                const register double tmp1_1 = C_1*w45;
1963                                const register double tmp2_1 = C_1*w48;
1964                                const register double tmp3_1 = C_0*w51;
1965                                const register double tmp4_1 = C_0*w44;
1966                                const register double tmp5_1 = C_1*w49;
1967                                const register double tmp6_1 = C_1*w50;
1968                                const register double tmp7_1 = C_0*w46;
1969                                EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1970                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
1971                                EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
1972                                EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
1973                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1974                                EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
1975                                EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
1976                                EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
1977                                EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
1978                                EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
1979                                EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
1980                                EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
1981                                EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
1982                                EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
1983                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
1984                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
1985                            }
1986                        }
1987                        ///////////////
1988                        // process D //
1989                        ///////////////
1990                        if (!D.isEmpty()) {
1991                            add_EM_S=true;
1992                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
1993                            if (D.actsExpanded()) {
1994                                const register double D_0 = D_p[0];
1995                                const register double D_1 = D_p[1];
1996                                const register double D_2 = D_p[2];
1997                                const register double D_3 = D_p[3];
1998                                const register double tmp0_0 = D_0 + D_1;
1999                                const register double tmp1_0 = D_2 + D_3;
2000                                const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;
2001                                const register double tmp3_0 = D_1 + D_2;
2002                                const register double tmp4_0 = D_1 + D_3;
2003                                const register double tmp5_0 = D_0 + D_2;
2004                                const register double tmp6_0 = D_0 + D_3;
2005                                const register double tmp0_1 = tmp0_0*w52;
2006                                const register double tmp1_1 = tmp1_0*w53;
2007                                const register double tmp2_1 = tmp2_0*w54;
2008                                const register double tmp3_1 = tmp1_0*w52;
2009                                const register double tmp4_1 = tmp0_0*w53;
2010                                const register double tmp5_1 = tmp3_0*w54;
2011                                const register double tmp6_1 = D_0*w55;
2012                                const register double tmp7_1 = D_3*w56;
2013                                const register double tmp8_1 = D_3*w55;
2014                                const register double tmp9_1 = D_0*w56;
2015                                const register double tmp10_1 = tmp4_0*w52;
2016                                const register double tmp11_1 = tmp5_0*w53;
2017                                const register double tmp12_1 = tmp5_0*w52;
2018                                const register double tmp13_1 = tmp4_0*w53;
2019                                const register double tmp14_1 = tmp6_0*w54;
2020                                const register double tmp15_1 = D_2*w55;
2021                                const register double tmp16_1 = D_1*w56;
2022                                const register double tmp17_1 = D_1*w55;
2023                                const register double tmp18_1 = D_2*w56;
2024                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2025                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2026                                EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
2027                                EM_S[INDEX2(3,0,4)]+=tmp2_1;
2028                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2029                                EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
2030                                EM_S[INDEX2(2,1,4)]+=tmp2_1;
2031                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
2032                                EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
2033                                EM_S[INDEX2(1,2,4)]+=tmp2_1;
2034                                EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
2035                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
2036                                EM_S[INDEX2(0,3,4)]+=tmp2_1;
2037                                EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
2038                                EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2039                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2040                            } else { // constant data
2041                                const register double tmp0_1 = D_p[0]*w57;
2042                                const register double tmp1_1 = D_p[0]*w58;
2043                                const register double tmp2_1 = D_p[0]*w59;
2044                                EM_S[INDEX2(0,0,4)]+=tmp2_1;
2045                                EM_S[INDEX2(1,0,4)]+=tmp0_1;
2046                                EM_S[INDEX2(2,0,4)]+=tmp0_1;
2047                                EM_S[INDEX2(3,0,4)]+=tmp1_1;
2048                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
2049                                EM_S[INDEX2(1,1,4)]+=tmp2_1;
2050                                EM_S[INDEX2(2,1,4)]+=tmp1_1;
2051                                EM_S[INDEX2(3,1,4)]+=tmp0_1;
2052                                EM_S[INDEX2(0,2,4)]+=tmp0_1;
2053                                EM_S[INDEX2(1,2,4)]+=tmp1_1;
2054                                EM_S[INDEX2(2,2,4)]+=tmp2_1;
2055                                EM_S[INDEX2(3,2,4)]+=tmp0_1;
2056                                EM_S[INDEX2(0,3,4)]+=tmp1_1;
2057                                EM_S[INDEX2(1,3,4)]+=tmp0_1;
2058                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
2059                                EM_S[INDEX2(3,3,4)]+=tmp2_1;
2060                            }
2061                        }
2062                        ///////////////
2063                        // process X //
2064                        ///////////////
2065                        if (!X.isEmpty()) {
2066                            add_EM_F=true;
2067                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2068                            if (X.actsExpanded()) {
2069                                const register double X_0_0 = X_p[INDEX2(0,0,2)];
2070                                const register double X_1_0 = X_p[INDEX2(1,0,2)];
2071                                const register double X_0_1 = X_p[INDEX2(0,1,2)];
2072                                const register double X_1_1 = X_p[INDEX2(1,1,2)];
2073                                const register double X_0_2 = X_p[INDEX2(0,2,2)];
2074                                const register double X_1_2 = X_p[INDEX2(1,2,2)];
2075                                const register double X_0_3 = X_p[INDEX2(0,3,2)];
2076                                const register double X_1_3 = X_p[INDEX2(1,3,2)];
2077                                const register double tmp0_0 = X_0_2 + X_0_3;
2078                                const register double tmp1_0 = X_1_1 + X_1_3;
2079                                const register double tmp2_0 = X_1_0 + X_1_2;
2080                                const register double tmp3_0 = X_0_0 + X_0_1;
2081                                const register double tmp0_1 = tmp0_0*w63;
2082                                const register double tmp1_1 = tmp1_0*w62;
2083                                const register double tmp2_1 = tmp2_0*w61;
2084                                const register double tmp3_1 = tmp3_0*w60;
2085                                const register double tmp4_1 = tmp0_0*w65;
2086                                const register double tmp5_1 = tmp3_0*w64;
2087                                const register double tmp6_1 = tmp2_0*w62;
2088                                const register double tmp7_1 = tmp1_0*w61;
2089                                const register double tmp8_1 = tmp2_0*w66;
2090                                const register double tmp9_1 = tmp3_0*w63;
2091                                const register double tmp10_1 = tmp0_0*w60;
2092                                const register double tmp11_1 = tmp1_0*w67;
2093                                const register double tmp12_1 = tmp1_0*w66;
2094                                const register double tmp13_1 = tmp3_0*w65;
2095                                const register double tmp14_1 = tmp0_0*w64;
2096                                const register double tmp15_1 = tmp2_0*w67;
2097                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2098                                EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2099                                EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2100                                EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2101                            } else { // constant data
2102                                const register double tmp0_1 = X_p[1]*w69;
2103                                const register double tmp1_1 = X_p[0]*w68;
2104                                const register double tmp2_1 = X_p[0]*w70;
2105                                const register double tmp3_1 = X_p[1]*w71;
2106                                EM_F[0]+=tmp0_1 + tmp1_1;
2107                                EM_F[1]+=tmp0_1 + tmp2_1;
2108                                EM_F[2]+=tmp1_1 + tmp3_1;
2109                                EM_F[3]+=tmp2_1 + tmp3_1;
2110                            }
2111                        }
2112                        ///////////////
2113                        // process Y //
2114                        ///////////////
2115                        if (!Y.isEmpty()) {
2116                            add_EM_F=true;
2117                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2118                            if (Y.actsExpanded()) {
2119                                const register double Y_0 = Y_p[0];
2120                                const register double Y_1 = Y_p[1];
2121                                const register double Y_2 = Y_p[2];
2122                                const register double Y_3 = Y_p[3];
2123                                const register double tmp0_0 = Y_1 + Y_2;
2124                                const register double tmp1_0 = Y_0 + Y_3;
2125                                const register double tmp0_1 = Y_0*w72;
2126                                const register double tmp1_1 = tmp0_0*w73;
2127                                const register double tmp2_1 = Y_3*w74;
2128                                const register double tmp3_1 = Y_1*w72;
2129                                const register double tmp4_1 = tmp1_0*w73;
2130                                const register double tmp5_1 = Y_2*w74;
2131                                const register double tmp6_1 = Y_2*w72;
2132                                const register double tmp7_1 = Y_1*w74;
2133                                const register double tmp8_1 = Y_3*w72;
2134                                const register double tmp9_1 = Y_0*w74;
2135                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2136                                EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2137                                EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2138                                EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2139                            } else { // constant data
2140                                const register double tmp0_1 = Y_p[0]*w75;
2141                                EM_F[0]+=tmp0_1;
2142                                EM_F[1]+=tmp0_1;
2143                                EM_F[2]+=tmp0_1;
2144                                EM_F[3]+=tmp0_1;
2145                            }
2146                        }
2147    
2148                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2149                        const index_t firstNode=m_N0*k1+k0;
2150                        IndexVector rowIndex;
2151                        rowIndex.push_back(m_dofMap[firstNode]);
2152                        rowIndex.push_back(m_dofMap[firstNode+1]);
2153                        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2154                        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2155                        if (add_EM_F) {
2156                            double *F_p=rhs.getSampleDataRW(0);
2157                            for (index_t i=0; i<rowIndex.size(); i++) {
2158                                if (rowIndex[i]<getNumDOF()) {
2159                                    F_p[rowIndex[i]]+=EM_F[i];
2160                                }
2161                            }
2162                        }
2163                        if (add_EM_S) {
2164                            addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2165                        }
2166    
2167                    } // end k0 loop
2168                } // end k1 loop
2169            } // end of colouring
2170        } // end of parallel region
2171    }
2172    
2173    //protected
2174    void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2175            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2176            const escript::Data& C, const escript::Data& D,
2177            const escript::Data& X, const escript::Data& Y,
2178            const escript::Data& d, const escript::Data& y) const
2179    {
2180        const double h0 = m_l0/m_gNE0;
2181        const double h1 = m_l1/m_gNE1;
2182        const double w0 = -.25*h1/h0;
2183        const double w1 = .25;
2184        const double w2 = -.25;
2185        const double w3 = .25*h0/h1;
2186        const double w4 = -.25*h0/h1;
2187        const double w5 = .25*h1/h0;
2188        const double w6 = -.125*h1;
2189        const double w7 = -.125*h0;
2190        const double w8 = .125*h1;
2191        const double w9 = .125*h0;
2192        const double w10 = .0625*h0*h1;
2193        const double w11 = -.5*h1;
2194        const double w12 = -.5*h0;
2195        const double w13 = .5*h1;
2196        const double w14 = .5*h0;
2197        const double w15 = .25*h0*h1;
2198    
2199        rhs.requireWrite();
2200    #pragma omp parallel
2201        {
2202            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2203    #pragma omp for
2204                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2205                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2206                        bool add_EM_S=false;
2207                        bool add_EM_F=false;
2208                        vector<double> EM_S(4*4, 0);
2209                        vector<double> EM_F(4, 0);
2210                        const index_t e = k0 + m_NE0*k1;
2211                        ///////////////
2212                        // process A //
2213                        ///////////////
2214                        if (!A.isEmpty()) {
2215                            add_EM_S=true;
2216                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2217                            const register double A_00 = A_p[INDEX2(0,0,2)];
2218                            const register double A_10 = A_p[INDEX2(1,0,2)];
2219                            const register double A_01 = A_p[INDEX2(0,1,2)];
2220                            const register double A_11 = A_p[INDEX2(1,1,2)];
2221                            const register double tmp0_0 = A_01 + A_10;
2222                            const register double tmp0_1 = A_11*w3;
2223                            const register double tmp1_1 = A_00*w0;
2224                            const register double tmp2_1 = A_01*w1;
2225                            const register double tmp3_1 = A_10*w2;
2226                            const register double tmp4_1 = tmp0_0*w1;
2227                            const register double tmp5_1 = A_11*w4;
2228                            const register double tmp6_1 = A_00*w5;
2229                            const register double tmp7_1 = tmp0_0*w2;
2230                            const register double tmp8_1 = A_10*w1;
2231                            const register double tmp9_1 = A_01*w2;
2232                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2233                            EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2234                            EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2235                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2236                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2237                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2238                            EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2239                            EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2240                            EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2241                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2242                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2243                            EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2244                            EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2245                            EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2246                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2247                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2248                        }
2249                        ///////////////
2250                        // process B //
2251                        ///////////////
2252                        if (!B.isEmpty()) {
2253                            add_EM_S=true;
2254                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2255                            const register double tmp2_1 = B_p[0]*w8;
2256                            const register double tmp0_1 = B_p[0]*w6;
2257                            const register double tmp3_1 = B_p[1]*w9;
2258                            const register double tmp1_1 = B_p[1]*w7;
2259                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2260                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2261                            EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
2262                            EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;
2263                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2264                            EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;
2265                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;
2266                            EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;
2267                            EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;
2268                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
2269                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;
2270                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2271                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;
2272                            EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;
2273                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;
2274                            EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
2275                        }
2276                        ///////////////
2277                        // process C //
2278                        ///////////////
2279                        if (!C.isEmpty()) {
2280                            add_EM_S=true;
2281                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2282                            const register double tmp1_1 = C_p[1]*w7;
2283                            const register double tmp0_1 = C_p[0]*w8;
2284                            const register double tmp3_1 = C_p[0]*w6;
2285                            const register double tmp2_1 = C_p[1]*w9;
2286                            EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2287                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2288                            EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
2289                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
2290                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2291                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;
2292                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;
2293                            EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;
2294                            EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;
2295                            EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
2296                            EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
2297                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2298                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;
2299                            EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;
2300                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2301                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;
2302                        }
2303                        ///////////////
2304                        // process D //
2305                        ///////////////
2306                        if (!D.isEmpty()) {
2307                            add_EM_S=true;
2308                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2309                            const register double tmp0_1 = D_p[0]*w10;
2310                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
2311                            EM_S[INDEX2(1,0,4)]+=tmp0_1;
2312                            EM_S[INDEX2(2,0,4)]+=tmp0_1;
2313                            EM_S[INDEX2(3,0,4)]+=tmp0_1;
2314                            EM_S[INDEX2(0,1,4)]+=tmp0_1;
2315                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
2316                            EM_S[INDEX2(2,1,4)]+=tmp0_1;
2317                            EM_S[INDEX2(3,1,4)]+=tmp0_1;
2318                            EM_S[INDEX2(0,2,4)]+=tmp0_1;
2319                            EM_S[INDEX2(1,2,4)]+=tmp0_1;
2320                            EM_S[INDEX2(2,2,4)]+=tmp0_1;
2321                            EM_S[INDEX2(3,2,4)]+=tmp0_1;
2322                            EM_S[INDEX2(0,3,4)]+=tmp0_1;
2323                            EM_S[INDEX2(1,3,4)]+=tmp0_1;
2324                            EM_S[INDEX2(2,3,4)]+=tmp0_1;
2325                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
2326                        }
2327                        ///////////////
2328                        // process X //
2329                        ///////////////
2330                        if (!X.isEmpty()) {
2331                            add_EM_F=true;
2332                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2333                            const register double tmp0_1 = X_p[0]*w11;
2334                            const register double tmp2_1 = X_p[0]*w13;
2335                            const register double tmp1_1 = X_p[1]*w12;
2336                            const register double tmp3_1 = X_p[1]*w14;
2337                            EM_F[0]+=tmp0_1 + tmp1_1;
2338                            EM_F[1]+=tmp1_1 + tmp2_1;
2339                            EM_F[2]+=tmp0_1 + tmp3_1;
2340                            EM_F[3]+=tmp2_1 + tmp3_1;
2341                        }
2342                        ///////////////
2343                        // process Y //
2344                        ///////////////
2345                        if (!Y.isEmpty()) {
2346                            add_EM_F=true;
2347                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2348                            const register double tmp0_1 = Y_p[0]*w15;
2349                            EM_F[0]+=tmp0_1;
2350                            EM_F[1]+=tmp0_1;
2351                            EM_F[2]+=tmp0_1;
2352                            EM_F[3]+=tmp0_1;
2353                        }
2354    
2355                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2356                        const index_t firstNode=m_N0*k1+k0;
2357                        IndexVector rowIndex;
2358                        rowIndex.push_back(m_dofMap[firstNode]);
2359                        rowIndex.push_back(m_dofMap[firstNode+1]);
2360                        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2361                        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2362                        if (add_EM_F) {
2363                            double *F_p=rhs.getSampleDataRW(0);
2364                            for (index_t i=0; i<rowIndex.size(); i++) {
2365                                if (rowIndex[i]<getNumDOF()) {
2366                                    F_p[rowIndex[i]]+=EM_F[i];
2367                                }
2368                            }
2369                        }
2370                        if (add_EM_S) {
2371                            addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2372                        }
2373    
2374                    } // end k0 loop
2375                } // end k1 loop
2376            } // end of colouring
2377        } // end of parallel region
2378    }
2379    
2380    //protected
2381    void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2382            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2383            const escript::Data& C, const escript::Data& D,
2384            const escript::Data& X, const escript::Data& Y,
2385            const escript::Data& d, const escript::Data& y) const
2386    {
2387        const double h0 = m_l0/m_gNE0;
2388        const double h1 = m_l1/m_gNE1;
2389        dim_t numEq, numComp;
2390        if (!mat)
2391            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
2392        else {
2393            numEq=mat->logical_row_block_size;
2394            numComp=mat->logical_col_block_size;
2395        }
2396    
2397        const double w0 = -0.1555021169820365539*h1/h0;
2398        const double w1 = 0.041666666666666666667;
2399        const double w2 = -0.15550211698203655390;
2400        const double w3 = 0.041666666666666666667*h0/h1;
2401        const double w4 = 0.15550211698203655390;
2402        const double w5 = -0.041666666666666666667;
2403        const double w6 = -0.01116454968463011277*h1/h0;
2404        const double w7 = 0.011164549684630112770;
2405        const double w8 = -0.011164549684630112770;
2406        const double w9 = -0.041666666666666666667*h1/h0;
2407        const double w10 = -0.041666666666666666667*h0/h1;
2408        const double w11 = 0.1555021169820365539*h1/h0;
2409        const double w12 = 0.1555021169820365539*h0/h1;
2410        const double w13 = 0.01116454968463011277*h0/h1;
2411        const double w14 = 0.01116454968463011277*h1/h0;
2412        const double w15 = 0.041666666666666666667*h1/h0;
2413        const double w16 = -0.01116454968463011277*h0/h1;
2414        const double w17 = -0.1555021169820365539*h0/h1;
2415        const double w18 = -0.33333333333333333333*h1/h0;
2416        const double w19 = 0.25000000000000000000;
2417        const double w20 = -0.25000000000000000000;
2418        const double w21 = 0.16666666666666666667*h0/h1;
2419        const double w22 = -0.16666666666666666667*h1/h0;
2420        const double w23 = -0.16666666666666666667*h0/h1;
2421        const double w24 = 0.33333333333333333333*h1/h0;
2422        const double w25 = 0.33333333333333333333*h0/h1;
2423        const double w26 = 0.16666666666666666667*h1/h0;
2424        const double w27 = -0.33333333333333333333*h0/h1;
2425        const double w28 = -0.032861463941450536761*h1;
2426        const double w29 = -0.032861463941450536761*h0;
2427        const double w30 = -0.12264065304058601714*h1;
2428        const double w31 = -0.0023593469594139828636*h1;
2429        const double w32 = -0.008805202725216129906*h0;
2430        const double w33 = -0.008805202725216129906*h1;
2431        const double w34 = 0.032861463941450536761*h1;
2432        const double w35 = 0.008805202725216129906*h1;
2433        const double w36 = 0.008805202725216129906*h0;
2434        const double w37 = 0.0023593469594139828636*h1;
2435        const double w38 = 0.12264065304058601714*h1;
2436        const double w39 = 0.032861463941450536761*h0;
2437        const double w40 = -0.12264065304058601714*h0;
2438        const double w41 = -0.0023593469594139828636*h0;
2439        const double w42 = 0.0023593469594139828636*h0;
2440        const double w43 = 0.12264065304058601714*h0;
2441        const double w44 = -0.16666666666666666667*h1;
2442        const double w45 = -0.083333333333333333333*h0;
2443        const double w46 = 0.083333333333333333333*h1;
2444        const double w47 = 0.16666666666666666667*h1;
2445        const double w48 = 0.083333333333333333333*h0;
2446        const double w49 = -0.16666666666666666667*h0;
2447        const double w50 = 0.16666666666666666667*h0;
2448        const double w51 = -0.083333333333333333333*h1;
2449        const double w52 = 0.025917019497006092316*h0*h1;
2450        const double w53 = 0.0018607582807716854616*h0*h1;
2451        const double w54 = 0.0069444444444444444444*h0*h1;
2452        const double w55 = 0.09672363354357992482*h0*h1;
2453        const double w56 = 0.00049858867864229740201*h0*h1;
2454        const double w57 = 0.055555555555555555556*h0*h1;
2455        const double w58 = 0.027777777777777777778*h0*h1;
2456        const double w59 = 0.11111111111111111111*h0*h1;
2457        const double w60 = -0.19716878364870322056*h1;
2458        const double w61 = -0.19716878364870322056*h0;
2459        const double w62 = -0.052831216351296779436*h0;
2460        const double w63 = -0.052831216351296779436*h1;
2461        const double w64 = 0.19716878364870322056*h1;
2462        const double w65 = 0.052831216351296779436*h1;
2463        const double w66 = 0.19716878364870322056*h0;
2464        const double w67 = 0.052831216351296779436*h0;
2465        const double w68 = -0.5*h1;
2466        const double w69 = -0.5*h0;
2467        const double w70 = 0.5*h1;
2468        const double w71 = 0.5*h0;
2469        const double w72 = 0.1555021169820365539*h0*h1;
2470        const double w73 = 0.041666666666666666667*h0*h1;
2471        const double w74 = 0.01116454968463011277*h0*h1;
2472        const double w75 = 0.25*h0*h1;
2473    
2474        rhs.requireWrite();
2475    #pragma omp parallel
2476        {
2477            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2478    #pragma omp for
2479                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2480                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2481                        bool add_EM_S=false;
2482                        bool add_EM_F=false;
2483                        vector<double> EM_S(4*4*numEq*numComp, 0);
2484                        vector<double> EM_F(4*numEq, 0);
2485                        const index_t e = k0 + m_NE0*k1;
2486                        ///////////////
2487                        // process A //
2488                        ///////////////
2489                        if (!A.isEmpty()) {
2490                            add_EM_S=true;
2491                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2492                            if (A.actsExpanded()) {
2493                                for (index_t k=0; k<numEq; k++) {
2494                                    for (index_t m=0; m<numComp; m++) {
2495                                        const register double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];
2496                                        const register double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];
2497                                        const register double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];
2498                                        const register double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];
2499                                        const register double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];
2500                                        const register double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];
2501                                        const register double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];
2502                                        const register double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];
2503                                        const register double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];
2504                                        const register double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];
2505                                        const register double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];
2506                                        const register double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];
2507                                        const register double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];
2508                                        const register double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];
2509                                        const register double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];
2510                                        const register double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];
2511                                        const register double tmp0_0 = A_01_0 + A_01_3;
2512                                        const register double tmp1_0 = A_00_0 + A_00_1;
2513                                        const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
2514                                        const register double tmp3_0 = A_00_2 + A_00_3;
2515                                        const register double tmp4_0 = A_10_1 + A_10_2;
2516                                        const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
2517                                        const register double tmp6_0 = A_01_3 + A_10_0;
2518                                        const register double tmp7_0 = A_01_0 + A_10_3;
2519                                        const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
2520                                        const register double tmp9_0 = A_01_0 + A_10_0;
2521                                        const register double tmp10_0 = A_01_3 + A_10_3;
2522                                        const register double tmp11_0 = A_11_1 + A_11_3;
2523                                        const register double tmp12_0 = A_11_0 + A_11_2;
2524                                        const register double tmp13_0 = A_01_2 + A_10_1;
2525                                        const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
2526                                        const register double tmp15_0 = A_01_1 + A_10_2;
2527                                        const register double tmp16_0 = A_10_0 + A_10_3;
2528                                        const register double tmp17_0 = A_01_1 + A_01_2;
2529                                        const register double tmp18_0 = A_01_1 + A_10_1;
2530                                        const register double tmp19_0 = A_01_2 + A_10_2;
2531                                        const register double tmp0_1 = A_10_3*w8;
2532                                        const register double tmp1_1 = tmp0_0*w1;
2533                                        const register double tmp2_1 = A_01_1*w4;
2534                                        const register double tmp3_1 = tmp1_0*w0;
2535                                        const register double tmp4_1 = A_01_2*w7;
2536                                        const register double tmp5_1 = tmp2_0*w3;
2537                                        const register double tmp6_1 = tmp3_0*w6;
2538                                        const register double tmp7_1 = A_10_0*w2;
2539                                        const register double tmp8_1 = tmp4_0*w5;
2540                                        const register double tmp9_1 = tmp2_0*w10;
2541                                        const register double tmp10_1 = tmp5_0*w9;
2542                                        const register double tmp11_1 = tmp6_0*w7;
2543                                        const register double tmp12_1 = tmp7_0*w4;
2544                                        const register double tmp13_1 = tmp8_0*w1;
2545                                        const register double tmp14_1 = A_10_0*w8;
2546                                        const register double tmp15_1 = A_01_2*w4;
2547                                        const register double tmp16_1 = tmp3_0*w0;
2548                                        const register double tmp17_1 = A_01_1*w7;
2549                                        const register double tmp18_1 = tmp1_0*w6;
2550                                        const register double tmp19_1 = A_10_3*w2;
2551                                        const register double tmp20_1 = tmp9_0*w4;
2552                                        const register double tmp21_1 = tmp1_0*w11;
2553                                        const register double tmp22_1 = tmp10_0*w7;
2554                                        const register double tmp23_1 = tmp3_0*w14;
2555                                        const register double tmp24_1 = tmp11_0*w13;
2556                                        const register double tmp25_1 = tmp12_0*w12;
2557                                        const register double tmp26_1 = tmp10_0*w4;
2558                                        const register double tmp27_1 = tmp3_0*w11;
2559                                        const register double tmp28_1 = tmp9_0*w7;
2560                                        const register double tmp29_1 = tmp1_0*w14;
2561                                        const register double tmp30_1 = tmp12_0*w13;
2562                                        const register double tmp31_1 = tmp11_0*w12;
2563                                        const register double tmp32_1 = tmp13_0*w2;
2564                                        const register double tmp33_1 = tmp14_0*w5;
2565                                        const register double tmp34_1 = tmp15_0*w8;
2566                                        const register double tmp35_1 = A_01_0*w8;
2567                                        const register double tmp36_1 = tmp16_0*w1;
2568                                        const register double tmp37_1 = A_10_1*w4;
2569                                        const register double tmp38_1 = tmp5_0*w15;
2570                                        const register double tmp39_1 = A_10_2*w7;
2571                                        const register double tmp40_1 = tmp11_0*w17;
2572                                        const register double tmp41_1 = A_01_3*w2;
2573                                        const register double tmp42_1 = tmp12_0*w16;
2574                                        const register double tmp43_1 = tmp17_0*w5;
2575                                        const register double tmp44_1 = tmp7_0*w7;
2576                                        const register double tmp45_1 = tmp6_0*w4;
2577                                        const register double tmp46_1 = A_01_3*w8;
2578                                        const register double tmp47_1 = A_10_2*w4;
2579                                        const register double tmp48_1 = A_10_1*w7;
2580                                        const register double tmp49_1 = tmp12_0*w17;
2581                                        const register double tmp50_1 = A_01_0*w2;
2582                                        const register double tmp51_1 = tmp11_0*w16;
2583                                        const register double tmp52_1 = tmp18_0*w8;
2584                                        const register double tmp53_1 = tmp19_0*w2;
2585                                        const register double tmp54_1 = tmp13_0*w8;
2586                                        const register double tmp55_1 = tmp15_0*w2;
2587                                        const register double tmp56_1 = tmp19_0*w8;
2588                                        const register double tmp57_1 = tmp18_0*w2;
2589                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
2590                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
2591                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
2592                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
2593                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;