/[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 3759 by caltinay, Fri Jan 6 06:54:51 2012 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 49  Rectangle::Rectangle(int n0, int n1, dou Line 49  Rectangle::Rectangle(int n0, int n1, dou
49      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
50          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
51    
52      // local number of elements (including overlap)      // local number of elements (with and without overlap)
53      m_NE0 = (m_NX>1 ? (n0+1)/m_NX : n0);      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)      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
55          m_NE0++;          m_NE0++;
56      m_NE1 = (m_NY>1 ? (n1+1)/m_NY : n1);      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)      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)
61          m_NE1++;          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      // local number of nodes
66      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
# Line 71  Rectangle::Rectangle(int n0, int n1, dou Line 76  Rectangle::Rectangle(int n0, int n1, dou
76    
77      populateSampleIds();      populateSampleIds();
78      createPattern();      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 106  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 126  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 141  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 241  const int* Rectangle::borrowSampleRefere Line 252  const int* Rectangle::borrowSampleRefere
252  {  {
253      switch (fsType) {      switch (fsType) {
254          case Nodes:          case Nodes:
         case ReducedNodes: //FIXME: reduced  
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:          case DegreesOfFreedom:
         case ReducedDegreesOfFreedom: //FIXME: reduced  
261              return &m_dofId[0];              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:          case ReducedElements:
268              return &m_elementId[0];              return &m_elementId[0];
# Line 269  bool Rectangle::ownSample(int fsType, in Line 286  bool Rectangle::ownSample(int fsType, in
286    
287      switch (fsType) {      switch (fsType) {
288          case Nodes:          case Nodes:
         case ReducedNodes: //FIXME: reduced  
289              return (m_dofMap[id] < getNumDOF());              return (m_dofMap[id] < getNumDOF());
290            case ReducedNodes:
291                if (m_coarseMesh)
292                    return m_coarseMesh->ownSample(Nodes, id);
293                break;
294          case DegreesOfFreedom:          case DegreesOfFreedom:
295          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
296              return true;              return true;
# Line 281  bool Rectangle::ownSample(int fsType, in Line 301  bool Rectangle::ownSample(int fsType, in
301          case FaceElements:          case FaceElements:
302          case ReducedFaceElements:          case ReducedFaceElements:
303              {              {
304                  // check ownership of face element's first node                  // determine which face the sample belongs to before
305                    // checking ownership of face element's first node
306                  const IndexVector faces = getNumFacesPerBoundary();                  const IndexVector faces = getNumFacesPerBoundary();
307                  dim_t n=0;                  dim_t n=0;
308                  for (size_t i=0; i<faces.size(); i++) {                  for (size_t i=0; i<faces.size(); i++) {
# Line 311  bool Rectangle::ownSample(int fsType, in Line 332  bool Rectangle::ownSample(int fsType, in
332      throw RipleyException(msg.str());      throw RipleyException(msg.str());
333  }  }
334    
335  void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const  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 {
431            stringstream msg;
432            msg << "setToNormal() not implemented for "
433                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
434            throw RipleyException(msg.str());
435        }
436    }
437    
438    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;
# Line 335  void Rectangle::setToGradient(escript::D Line 639  void Rectangle::setToGradient(escript::D
639      const double cy7 = 1./h1;      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();
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 356  void Rectangle::setToGradient(escript::D Line 660  void Rectangle::setToGradient(escript::D
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();
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 373  void Rectangle::setToGradient(escript::D Line 676  void Rectangle::setToGradient(escript::D
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) {
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) {
# Line 442  void Rectangle::setToGradient(escript::D Line 744  void Rectangle::setToGradient(escript::D
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) {
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) {
# Line 504  void Rectangle::setToGradient(escript::D Line 805  void Rectangle::setToGradient(escript::D
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 544  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 563  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 626  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 636  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 646  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 656  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 668  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());  
     }  
 }  
   
 void Rectangle::setToNormal(escript::Data& out) const  
 {  
     if (out.getFunctionSpace().getTypeCode() == FaceElements) {  
 #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);  
                     // set vector at two quadrature points  
                     *o++ = -1.;  
                     *o++ = 0.;  
                     *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);  
                     // 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.;  
                 }  
             }  
   
             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());  
     }  
 }  
   
 Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  
                                                 bool reducedColOrder) const  
 {  
     if (reducedRowOrder || reducedColOrder)  
         throw RipleyException("getPattern() not implemented for reduced order");  
   
     return m_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;  
         }  
     }  
 }  
   
 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::getNumDOF() const  
 {  
     return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;  
 }  
   
 //protected  
 dim_t Rectangle::getNumFaceElements() const  
 {  
     const IndexVector faces = getNumFacesPerBoundary();  
     dim_t n=0;  
     for (size_t i=0; i<faces.size(); i++)  
         n+=faces[i];  
     return n;  
 }  
   
 //protected  
 void Rectangle::assembleCoordinates(escript::Data& arg) const  
 {  
     escriptDataC x = arg.getDataC();  
     int numDim = m_numDim;  
     if (!isDataPointShapeEqual(&x, 1, &numDim))  
         throw RipleyException("setToX: Invalid Data object shape");  
     if (!numSamplesEqual(&x, 1, getNumNodes()))  
         throw RipleyException("setToX: Illegal number of samples in Data object");  
   
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
     arg.requireWrite();  
 #pragma omp parallel for  
     for (dim_t i1 = 0; i1 < m_N1; i1++) {  
         for (dim_t i0 = 0; i0 < m_N0; i0++) {  
             double* point = arg.getSampleDataRW(i0+m_N0*i1);  
             point[0] = xdx.first+i0*xdx.second;  
             point[1] = ydy.first+i1*ydy.second;  
         }  
968      }      }
969  }  }
970    
# Line 1041  void Rectangle::createPattern() Line 1123  void Rectangle::createPattern()
1123      // The rest is assigned in the loop further down      // The rest is assigned in the loop further down
1124      m_dofMap.assign(getNumNodes(), 0);      m_dofMap.assign(getNumNodes(), 0);
1125  #pragma omp parallel for  #pragma omp parallel for
1126      for (index_t i=bottom; i<m_N1; i++) {      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1127          for (index_t j=left; j<m_N0; j++) {          for (index_t j=left; j<left+nDOF0; j++) {
1128              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
1129          }          }
1130      }      }
# Line 1209  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 c0 = .25;          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) {
# Line 1224  void Rectangle::interpolateNodesOnElemen Line 1306  void Rectangle::interpolateNodesOnElemen
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 c0 = .16666666666666666667;          const double c0 = .16666666666666666667;
1312          const double c1 = .044658198738520451079;          const double c1 = .044658198738520451079;
1313          const double c2 = .62200846792814621559;          const double c2 = .62200846792814621559;
# Line 1246  void Rectangle::interpolateNodesOnElemen Line 1327  void Rectangle::interpolateNodesOnElemen
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 1256  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            out.requireWrite();
1340          const double c0 = .5;          const double c0 = .5;
1341  #pragma omp parallel  #pragma omp parallel
1342          {          {
             /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */  
1343              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1344  #pragma omp for nowait  #pragma omp for nowait
1345                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 1304  void Rectangle::interpolateNodesOnFaces( Line 1384  void Rectangle::interpolateNodesOnFaces(
1384                      } /* end of component loop i */                      } /* end of component loop i */
1385                  } /* end of k0 loop */                  } /* end of k0 loop */
1386              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */  
1387          } // end of parallel section          } // end of parallel section
1388      } else {      } else {
1389            out.requireWrite();
1390          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1391          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1392  #pragma omp parallel  #pragma omp parallel
1393          {          {
             /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */  
1394              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1395  #pragma omp for nowait  #pragma omp for nowait
1396                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 1360  void Rectangle::interpolateNodesOnFaces( Line 1439  void Rectangle::interpolateNodesOnFaces(
1439                      } /* end of component loop i */                      } /* end of component loop i */
1440                  } /* end of k0 loop */                  } /* end of k0 loop */
1441              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */  
1442          } // end of parallel section          } // end of parallel section
1443      }      }
1444  }  }
# Line 1374  void Rectangle::assemblePDESingle(Paso_S Line 1452  void Rectangle::assemblePDESingle(Paso_S
1452  {  {
1453      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1454      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
     /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */  
1455      const double w0 = -0.1555021169820365539*h1/h0;      const double w0 = -0.1555021169820365539*h1/h0;
1456      const double w1 = 0.041666666666666666667;      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;      const double w10 = -0.041666666666666666667*h0/h1;
1466      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = 0.1555021169820365539*h1/h0;
1467      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = 0.1555021169820365539*h0/h1;
# Line 1387  void Rectangle::assemblePDESingle(Paso_S Line 1472  void Rectangle::assemblePDESingle(Paso_S
1472      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*h0/h1;
1473      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*h1/h0;
1474      const double w19 = 0.25000000000000000000;      const double w19 = 0.25000000000000000000;
     const double w2 = -0.15550211698203655390;  
1475      const double w20 = -0.25000000000000000000;      const double w20 = -0.25000000000000000000;
1476      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*h0/h1;
1477      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*h1/h0;
# Line 1398  void Rectangle::assemblePDESingle(Paso_S Line 1482  void Rectangle::assemblePDESingle(Paso_S
1482      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = -0.33333333333333333333*h0/h1;
1483      const double w28 = -0.032861463941450536761*h1;      const double w28 = -0.032861463941450536761*h1;
1484      const double w29 = -0.032861463941450536761*h0;      const double w29 = -0.032861463941450536761*h0;
     const double w3 = 0.041666666666666666667*h0/h1;  
1485      const double w30 = -0.12264065304058601714*h1;      const double w30 = -0.12264065304058601714*h1;
1486      const double w31 = -0.0023593469594139828636*h1;      const double w31 = -0.0023593469594139828636*h1;
1487      const double w32 = -0.008805202725216129906*h0;      const double w32 = -0.008805202725216129906*h0;
# Line 1409  void Rectangle::assemblePDESingle(Paso_S Line 1492  void Rectangle::assemblePDESingle(Paso_S
1492      const double w37 = 0.0023593469594139828636*h1;      const double w37 = 0.0023593469594139828636*h1;
1493      const double w38 = 0.12264065304058601714*h1;      const double w38 = 0.12264065304058601714*h1;
1494      const double w39 = 0.032861463941450536761*h0;      const double w39 = 0.032861463941450536761*h0;
     const double w4 = 0.15550211698203655390;  
1495      const double w40 = -0.12264065304058601714*h0;      const double w40 = -0.12264065304058601714*h0;
1496      const double w41 = -0.0023593469594139828636*h0;      const double w41 = -0.0023593469594139828636*h0;
1497      const double w42 = 0.0023593469594139828636*h0;      const double w42 = 0.0023593469594139828636*h0;
# Line 1420  void Rectangle::assemblePDESingle(Paso_S Line 1502  void Rectangle::assemblePDESingle(Paso_S
1502      const double w47 = 0.16666666666666666667*h1;      const double w47 = 0.16666666666666666667*h1;
1503      const double w48 = 0.083333333333333333333*h0;      const double w48 = 0.083333333333333333333*h0;
1504      const double w49 = -0.16666666666666666667*h0;      const double w49 = -0.16666666666666666667*h0;
     const double w5 = -0.041666666666666666667;  
1505      const double w50 = 0.16666666666666666667*h0;      const double w50 = 0.16666666666666666667*h0;
1506      const double w51 = -0.083333333333333333333*h1;      const double w51 = -0.083333333333333333333*h1;
1507      const double w52 = 0.025917019497006092316*h0*h1;      const double w52 = 0.025917019497006092316*h0*h1;
# Line 1431  void Rectangle::assemblePDESingle(Paso_S Line 1512  void Rectangle::assemblePDESingle(Paso_S
1512      const double w57 = 0.055555555555555555556*h0*h1;      const double w57 = 0.055555555555555555556*h0*h1;
1513      const double w58 = 0.027777777777777777778*h0*h1;      const double w58 = 0.027777777777777777778*h0*h1;
1514      const double w59 = 0.11111111111111111111*h0*h1;      const double w59 = 0.11111111111111111111*h0*h1;
     const double w6 = -0.01116454968463011277*h1/h0;  
1515      const double w60 = -0.19716878364870322056*h1;      const double w60 = -0.19716878364870322056*h1;
1516      const double w61 = -0.19716878364870322056*h0;      const double w61 = -0.19716878364870322056*h0;
1517      const double w62 = -0.052831216351296779436*h0;      const double w62 = -0.052831216351296779436*h0;
# Line 1442  void Rectangle::assemblePDESingle(Paso_S Line 1522  void Rectangle::assemblePDESingle(Paso_S
1522      const double w67 = 0.052831216351296779436*h0;      const double w67 = 0.052831216351296779436*h0;
1523      const double w68 = -0.5*h1;      const double w68 = -0.5*h1;
1524      const double w69 = -0.5*h0;      const double w69 = -0.5*h0;
     const double w7 = 0.011164549684630112770;  
1525      const double w70 = 0.5*h1;      const double w70 = 0.5*h1;
1526      const double w71 = 0.5*h0;      const double w71 = 0.5*h0;
1527      const double w72 = 0.1555021169820365539*h0*h1;      const double w72 = 0.1555021169820365539*h0*h1;
1528      const double w73 = 0.041666666666666666667*h0*h1;      const double w73 = 0.041666666666666666667*h0*h1;
1529      const double w74 = 0.01116454968463011277*h0*h1;      const double w74 = 0.01116454968463011277*h0*h1;
1530      const double w75 = 0.25*h0*h1;      const double w75 = 0.25*h0*h1;
     const double w8 = -0.011164549684630112770;  
     const double w9 = -0.041666666666666666667*h1/h0;  
     /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */  
1531    
1532      rhs.requireWrite();      rhs.requireWrite();
1533  #pragma omp parallel  #pragma omp parallel
# Line 1465  void Rectangle::assemblePDESingle(Paso_S Line 1541  void Rectangle::assemblePDESingle(Paso_S
1541                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1542                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1543                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
                     /*** GENERATOR SNIP_PDE_SINGLE TOP */  
1544                      ///////////////                      ///////////////
1545                      // process A //                      // process A //
1546                      ///////////////                      ///////////////
# Line 1474  void Rectangle::assemblePDESingle(Paso_S Line 1549  void Rectangle::assemblePDESingle(Paso_S
1549                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1550                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
1551                              const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];                              const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
                             const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];  
1552                              const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];                              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)];                              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)];                              const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
                             const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];  
1556                              const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];                              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)];                              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)];                              const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
                             const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];  
1560                              const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];                              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)];                              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)];                              const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
                             const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];  
1564                              const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];                              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)];                              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;                              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;                              const register double tmp12_0 = A_11_0 + A_11_2;
                             const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;  
1578                              const register double tmp10_0 = A_01_3 + A_10_3;                              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;                              const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
                             const register double tmp0_0 = A_01_0 + A_01_3;  
1580                              const register double tmp13_0 = A_01_2 + A_10_1;                              const register double tmp13_0 = A_01_2 + A_10_1;
                             const register double tmp3_0 = A_00_2 + A_00_3;  
1581                              const register double tmp11_0 = A_11_1 + A_11_3;                              const register double tmp11_0 = A_11_1 + A_11_3;
1582                              const register double tmp18_0 = A_01_1 + A_10_1;                              const register double tmp18_0 = A_01_1 + A_10_1;
                             const register double tmp1_0 = A_00_0 + A_00_1;  
1583                              const register double tmp15_0 = A_01_1 + A_10_2;                              const register double tmp15_0 = A_01_1 + A_10_2;
                             const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;  
1584                              const register double tmp16_0 = A_10_0 + A_10_3;                              const register double tmp16_0 = A_10_0 + A_10_3;
                             const register double tmp6_0 = A_01_3 + A_10_0;  
1585                              const register double tmp17_0 = A_01_1 + A_01_2;                              const register double tmp17_0 = A_01_1 + A_01_2;
                             const register double tmp9_0 = A_01_0 + A_10_0;  
                             const register double tmp7_0 = A_01_0 + A_10_3;  
                             const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;  
1586                              const register double tmp19_0 = A_01_2 + A_10_2;                              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;                              const register double tmp14_1 = A_10_0*w8;
1598                              const register double tmp23_1 = tmp3_0*w14;                              const register double tmp23_1 = tmp3_0*w14;
1599                              const register double tmp35_1 = A_01_0*w8;                              const register double tmp35_1 = A_01_0*w8;
1600                              const register double tmp54_1 = tmp13_0*w8;                              const register double tmp54_1 = tmp13_0*w8;
1601                              const register double tmp20_1 = tmp9_0*w4;                              const register double tmp20_1 = tmp9_0*w4;
1602                              const register double tmp25_1 = tmp12_0*w12;                              const register double tmp25_1 = tmp12_0*w12;
                             const register double tmp2_1 = A_01_1*w4;  
1603                              const register double tmp44_1 = tmp7_0*w7;                              const register double tmp44_1 = tmp7_0*w7;
1604                              const register double tmp26_1 = tmp10_0*w4;                              const register double tmp26_1 = tmp10_0*w4;
1605                              const register double tmp52_1 = tmp18_0*w8;                              const register double tmp52_1 = tmp18_0*w8;
1606                              const register double tmp48_1 = A_10_1*w7;                              const register double tmp48_1 = A_10_1*w7;
1607                              const register double tmp46_1 = A_01_3*w8;                              const register double tmp46_1 = A_01_3*w8;
1608                              const register double tmp50_1 = A_01_0*w2;                              const register double tmp50_1 = A_01_0*w2;
                             const register double tmp8_1 = tmp4_0*w5;  
1609                              const register double tmp56_1 = tmp19_0*w8;                              const register double tmp56_1 = tmp19_0*w8;
                             const register double tmp9_1 = tmp2_0*w10;  
1610                              const register double tmp19_1 = A_10_3*w2;                              const register double tmp19_1 = A_10_3*w2;
1611                              const register double tmp47_1 = A_10_2*w4;                              const register double tmp47_1 = A_10_2*w4;
1612                              const register double tmp16_1 = tmp3_0*w0;                              const register double tmp16_1 = tmp3_0*w0;
# Line 1537  void Rectangle::assemblePDESingle(Paso_S Line 1619  void Rectangle::assemblePDESingle(Paso_S
1619                              const register double tmp34_1 = tmp15_0*w8;                              const register double tmp34_1 = tmp15_0*w8;
1620                              const register double tmp33_1 = tmp14_0*w5;                              const register double tmp33_1 = tmp14_0*w5;
1621                              const register double tmp24_1 = tmp11_0*w13;                              const register double tmp24_1 = tmp11_0*w13;
                             const register double tmp3_1 = tmp1_0*w0;  
                             const register double tmp5_1 = tmp2_0*w3;  
1622                              const register double tmp43_1 = tmp17_0*w5;                              const register double tmp43_1 = tmp17_0*w5;
1623                              const register double tmp15_1 = A_01_2*w4;                              const register double tmp15_1 = A_01_2*w4;
1624                              const register double tmp53_1 = tmp19_0*w2;                              const register double tmp53_1 = tmp19_0*w2;
# Line 1557  void Rectangle::assemblePDESingle(Paso_S Line 1637  void Rectangle::assemblePDESingle(Paso_S
1637                              const register double tmp42_1 = tmp12_0*w16;                              const register double tmp42_1 = tmp12_0*w16;
1638                              const register double tmp49_1 = tmp12_0*w17;                              const register double tmp49_1 = tmp12_0*w17;
1639                              const register double tmp21_1 = tmp1_0*w11;                              const register double tmp21_1 = tmp1_0*w11;
                             const register double tmp1_1 = tmp0_0*w1;  
1640                              const register double tmp45_1 = tmp6_0*w4;                              const register double tmp45_1 = tmp6_0*w4;
                             const register double tmp7_1 = A_10_0*w2;  
                             const register double tmp6_1 = tmp3_0*w6;  
1641                              const register double tmp13_1 = tmp8_0*w1;                              const register double tmp13_1 = tmp8_0*w1;
1642                              const register double tmp36_1 = tmp16_0*w1;                              const register double tmp36_1 = tmp16_0*w1;
1643                              const register double tmp41_1 = A_01_3*w2;                              const register double tmp41_1 = A_01_3*w2;
1644                              const register double tmp30_1 = tmp12_0*w13;                              const register double tmp30_1 = tmp12_0*w13;
                             const register double tmp4_1 = A_01_2*w7;  
                             const register double tmp0_1 = A_10_3*w8;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
1645                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1646                              EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              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;                              EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
1649                              EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;                              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;                              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;                              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(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
                             EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
1655                              EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;                              EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
1656                              EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;                              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;                              EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
1658                              EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;                              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                          } else { /* constant data */                              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)];                              const register double A_00 = A_p[INDEX2(0,0,2)];
                             const register double A_01 = A_p[INDEX2(0,1,2)];  
1663                              const register double A_10 = A_p[INDEX2(1,0,2)];                              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)];                              const register double A_11 = A_p[INDEX2(1,1,2)];
1666                              const register double tmp0_0 = A_01 + A_10;                              const register double tmp0_0 = A_01 + A_10;
1667                              const register double tmp0_1 = A_00*w18;                              const register double tmp0_1 = A_00*w18;
1668                              const register double tmp10_1 = A_01*w20;                              const register double tmp1_1 = A_01*w19;
1669                              const register double tmp12_1 = A_00*w26;                              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;                              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;                              const register double tmp8_1 = A_00*w24;
                             const register double tmp13_1 = A_10*w19;  
1676                              const register double tmp9_1 = tmp0_0*w20;                              const register double tmp9_1 = tmp0_0*w20;
1677                              const register double tmp3_1 = A_11*w21;                              const register double tmp10_1 = A_01*w20;
1678                              const register double tmp11_1 = A_11*w27;                              const register double tmp11_1 = A_11*w27;
1679                              const register double tmp1_1 = A_01*w19;                              const register double tmp12_1 = A_00*w26;
1680                              const register double tmp6_1 = A_11*w23;                              const register double tmp13_1 = A_10*w19;
                             const register double tmp7_1 = A_11*w25;  
                             const register double tmp2_1 = A_10*w20;  
                             const register double tmp5_1 = tmp0_0*w19;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
1681                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1682                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              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;                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1685                              EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;                              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;                              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;                              EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1690                              EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
                             EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
1691                              EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1692                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                              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;                              EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1694                              EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;                              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                      ///////////////                      ///////////////
# Line 1636  void Rectangle::assemblePDESingle(Paso_S Line 1711  void Rectangle::assemblePDESingle(Paso_S
1711                              const register double B_1_2 = B_p[INDEX2(1,2,2)];                              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)];                              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)];                              const register double B_1_3 = B_p[INDEX2(1,3,2)];
1714                              const register double tmp3_0 = B_0_0 + B_0_2;                              const register double tmp0_0 = B_1_0 + B_1_1;
1715                              const register double tmp1_0 = B_1_2 + B_1_3;                              const register double tmp1_0 = B_1_2 + B_1_3;
1716                              const register double tmp2_0 = B_0_1 + B_0_3;                              const register double tmp2_0 = B_0_1 + B_0_3;
1717                              const register double tmp0_0 = B_1_0 + B_1_1;                              const register double tmp3_0 = B_0_0 + B_0_2;
1718                              const register double tmp63_1 = B_1_1*w42;                              const register double tmp63_1 = B_1_1*w42;
1719                              const register double tmp79_1 = B_1_1*w40;                              const register double tmp79_1 = B_1_1*w40;
1720                              const register double tmp37_1 = tmp3_0*w35;                              const register double tmp37_1 = tmp3_0*w35;
# Line 1720  void Rectangle::assemblePDESingle(Paso_S Line 1795  void Rectangle::assemblePDESingle(Paso_S
1795                              const register double tmp78_1 = B_0_0*w34;                              const register double tmp78_1 = B_0_0*w34;
1796                              const register double tmp12_1 = tmp0_0*w36;                              const register double tmp12_1 = tmp0_0*w36;
1797                              const register double tmp75_1 = B_1_0*w32;                              const register double tmp75_1 = B_1_0*w32;
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
1798                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              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(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              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;                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
1802                              EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;                              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;                              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;                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1807                              EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
1808                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;                              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(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                              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;                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
1811                              EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;                              EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1812                          } else { /* constant data */                              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];                              const register double B_0 = B_p[0];
1816                              const register double B_1 = B_p[1];                              const register double B_1 = B_p[1];
                             const register double tmp6_1 = B_1*w50;  
                             const register double tmp1_1 = B_1*w45;  
                             const register double tmp5_1 = B_1*w49;  
                             const register double tmp4_1 = B_1*w48;  
1817                              const register double tmp0_1 = B_0*w44;                              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;                              const register double tmp2_1 = B_0*w46;
                             const register double tmp7_1 = B_0*w51;  
1820                              const register double tmp3_1 = B_0*w47;                              const register double tmp3_1 = B_0*w47;
1821                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const register double tmp4_1 = B_1*w48;
1822                              EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;                              const register double tmp5_1 = B_1*w49;
1823                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;                              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;                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1826                              EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;                              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;                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
1829                              EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;                              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;                              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;                              EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
1834                              EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;                              EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
                             EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;  
1835                              EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;                              EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
1836                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1837                              EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;                              EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
1838                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;                              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                      ///////////////                      ///////////////
# Line 1780  void Rectangle::assemblePDESingle(Paso_S Line 1855  void Rectangle::assemblePDESingle(Paso_S
1855                              const register double C_1_2 = C_p[INDEX2(1,2,2)];                              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)];                              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)];                              const register double C_1_3 = C_p[INDEX2(1,3,2)];
1858                              const register double tmp2_0 = C_0_1 + C_0_3;                              const register double tmp0_0 = C_1_0 + C_1_1;
1859                              const register double tmp1_0 = C_1_2 + C_1_3;                              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;                              const register double tmp3_0 = C_0_0 + C_0_2;
                             const register double tmp0_0 = C_1_0 + C_1_1;  
1862                              const register double tmp64_1 = C_0_2*w30;                              const register double tmp64_1 = C_0_2*w30;
1863                              const register double tmp14_1 = C_0_2*w28;                              const register double tmp14_1 = C_0_2*w28;
1864                              const register double tmp19_1 = C_0_3*w31;                              const register double tmp19_1 = C_0_3*w31;
# Line 1864  void Rectangle::assemblePDESingle(Paso_S Line 1939  void Rectangle::assemblePDESingle(Paso_S
1939                              const register double tmp78_1 = C_0_0*w34;                              const register double tmp78_1 = C_0_0*w34;
1940                              const register double tmp11_1 = tmp0_0*w36;                              const register double tmp11_1 = tmp0_0*w36;
1941                              const register double tmp23_1 = C_0_2*w33;                              const register double tmp23_1 = C_0_2*w33;
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
1942                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              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(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              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;                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
1946                              EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;                              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;                              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;                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1951                              EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
1952                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;                              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(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                              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;                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
1955                              EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;                              EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1956                          } else { /* constant data */                              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];                              const register double C_0 = C_p[0];
1960                              const register double C_1 = C_p[1];                              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;                              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;                              const register double tmp3_1 = C_0*w51;
1965                              const register double tmp4_1 = C_0*w44;                              const register double tmp4_1 = C_0*w44;
                             const register double tmp7_1 = C_0*w46;  
1966                              const register double tmp5_1 = C_1*w49;                              const register double tmp5_1 = C_1*w49;
                             const register double tmp2_1 = C_1*w48;  
                             const register double tmp0_1 = C_0*w47;  
1967                              const register double tmp6_1 = C_1*w50;                              const register double tmp6_1 = C_1*w50;
1968                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const register double tmp7_1 = C_0*w46;
                             EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;  
1969                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1970                              EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;                              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;                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
1973                              EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;                              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;                              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;                              EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
1978                              EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
                             EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
1979                              EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;                              EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
1980                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;                              EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
1981                              EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;                              EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
1982                              EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;                              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                      ///////////////                      ///////////////
# Line 1920  void Rectangle::assemblePDESingle(Paso_S Line 1995  void Rectangle::assemblePDESingle(Paso_S
1995                              const register double D_1 = D_p[1];                              const register double D_1 = D_p[1];
1996                              const register double D_2 = D_p[2];                              const register double D_2 = D_p[2];
1997                              const register double D_3 = D_p[3];                              const register double D_3 = D_p[3];
                             const register double tmp4_0 = D_1 + D_3;  
                             const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;  
                             const register double tmp5_0 = D_0 + D_2;  
1998                              const register double tmp0_0 = D_0 + D_1;                              const register double tmp0_0 = D_0 + D_1;
                             const register double tmp6_0 = D_0 + D_3;  
1999                              const register double tmp1_0 = D_2 + D_3;                              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;                              const register double tmp3_0 = D_1 + D_2;
2002                              const register double tmp16_1 = D_1*w56;                              const register double tmp4_0 = D_1 + D_3;
2003                              const register double tmp14_1 = tmp6_0*w54;                              const register double tmp5_0 = D_0 + D_2;
2004                              const register double tmp8_1 = D_3*w55;                              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;                              const register double tmp2_1 = tmp2_0*w54;
                             const register double tmp12_1 = tmp5_0*w52;  
                             const register double tmp4_1 = tmp0_0*w53;  
2008                              const register double tmp3_1 = tmp1_0*w52;                              const register double tmp3_1 = tmp1_0*w52;
2009                              const register double tmp13_1 = tmp4_0*w53;                              const register double tmp4_1 = tmp0_0*w53;
                             const register double tmp10_1 = tmp4_0*w52;  
                             const register double tmp1_1 = tmp1_0*w53;  
                             const register double tmp7_1 = D_3*w56;  
                             const register double tmp0_1 = tmp0_0*w52;  
                             const register double tmp11_1 = tmp5_0*w53;  
                             const register double tmp9_1 = D_0*w56;  
2010                              const register double tmp5_1 = tmp3_0*w54;                              const register double tmp5_1 = tmp3_0*w54;
                             const register double tmp18_1 = D_2*w56;  
                             const register double tmp17_1 = D_1*w55;  
2011                              const register double tmp6_1 = D_0*w55;                              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;                              const register double tmp15_1 = D_2*w55;
2021                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const register double tmp16_1 = D_1*w56;
2022                              EM_S[INDEX2(1,2,4)]+=tmp2_1;                              const register double tmp17_1 = D_1*w55;
2023                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;                              const register double tmp18_1 = D_2*w56;
2024                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2025                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;                              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;                              EM_S[INDEX2(3,0,4)]+=tmp2_1;
2028                              EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;                              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;                              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;                              EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
2033                              EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;                              EM_S[INDEX2(1,2,4)]+=tmp2_1;
                             EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;  
2034                              EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;                              EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
2035                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
2036                              EM_S[INDEX2(0,3,4)]+=tmp2_1;                              EM_S[INDEX2(0,3,4)]+=tmp2_1;
2037                              EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;                              EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
2038                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2039                              const register double D_0 = D_p[0];                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2040                              const register double tmp2_1 = D_0*w59;                          } else { // constant data
2041                              const register double tmp1_1 = D_0*w58;                              const register double tmp0_1 = D_p[0]*w57;
2042                              const register double tmp0_1 = D_0*w57;                              const register double tmp1_1 = D_p[0]*w58;
2043                              EM_S[INDEX2(0,1,4)]+=tmp0_1;                              const register double tmp2_1 = D_p[0]*w59;
                             EM_S[INDEX2(1,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1;  
2044                              EM_S[INDEX2(0,0,4)]+=tmp2_1;                              EM_S[INDEX2(0,0,4)]+=tmp2_1;
2045                              EM_S[INDEX2(3,3,4)]+=tmp2_1;                              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;                              EM_S[INDEX2(3,0,4)]+=tmp1_1;
2048                              EM_S[INDEX2(3,1,4)]+=tmp0_1;                              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;                              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;                              EM_S[INDEX2(0,2,4)]+=tmp0_1;
2053                              EM_S[INDEX2(2,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,2,4)]+=tmp1_1;
                             EM_S[INDEX2(1,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1;  
2054                              EM_S[INDEX2(2,2,4)]+=tmp2_1;                              EM_S[INDEX2(2,2,4)]+=tmp2_1;
2055                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(3,2,4)]+=tmp0_1;
2056                              EM_S[INDEX2(0,3,4)]+=tmp1_1;                              EM_S[INDEX2(0,3,4)]+=tmp1_1;
2057                              EM_S[INDEX2(1,1,4)]+=tmp2_1;                              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                      ///////////////                      ///////////////
# Line 2000  void Rectangle::assemblePDESingle(Paso_S Line 2074  void Rectangle::assemblePDESingle(Paso_S
2074                              const register double X_1_2 = X_p[INDEX2(1,2,2)];                              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)];                              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)];                              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;                              const register double tmp1_0 = X_1_1 + X_1_3;
                             const register double tmp3_0 = X_0_0 + X_0_1;  
2079                              const register double tmp2_0 = X_1_0 + X_1_2;                              const register double tmp2_0 = X_1_0 + X_1_2;
2080                              const register double tmp0_0 = X_0_2 + X_0_3;                              const register double tmp3_0 = X_0_0 + X_0_1;
2081                              const register double tmp8_1 = tmp2_0*w66;                              const register double tmp0_1 = tmp0_0*w63;
2082                              const register double tmp5_1 = tmp3_0*w64;                              const register double tmp1_1 = tmp1_0*w62;
2083                              const register double tmp14_1 = tmp0_0*w64;                              const register double tmp2_1 = tmp2_0*w61;
2084                              const register double tmp3_1 = tmp3_0*w60;                              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;                              const register double tmp9_1 = tmp3_0*w63;
                             const register double tmp13_1 = tmp3_0*w65;  
                             const register double tmp12_1 = tmp1_0*w66;  
2091                              const register double tmp10_1 = tmp0_0*w60;                              const register double tmp10_1 = tmp0_0*w60;
                             const register double tmp2_1 = tmp2_0*w61;  
                             const register double tmp6_1 = tmp2_0*w62;  
                             const register double tmp4_1 = tmp0_0*w65;  
2092                              const register double tmp11_1 = tmp1_0*w67;                              const register double tmp11_1 = tmp1_0*w67;
2093                              const register double tmp1_1 = tmp1_0*w62;                              const register double tmp12_1 = tmp1_0*w66;
2094                              const register double tmp7_1 = tmp1_0*w61;                              const register double tmp13_1 = tmp3_0*w65;
2095                              const register double tmp0_1 = tmp0_0*w63;                              const register double tmp14_1 = tmp0_0*w64;
2096                              const register double tmp15_1 = tmp2_0*w67;                              const register double tmp15_1 = tmp2_0*w67;
2097                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2098                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2099                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2100                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2101                          } else { /* constant data */                          } else { // constant data
2102                              const register double X_0 = X_p[0];                              const register double tmp0_1 = X_p[1]*w69;
2103                              const register double X_1 = X_p[1];                              const register double tmp1_1 = X_p[0]*w68;
2104                              const register double tmp3_1 = X_1*w71;                              const register double tmp2_1 = X_p[0]*w70;
2105                              const register double tmp2_1 = X_0*w70;                              const register double tmp3_1 = X_p[1]*w71;
                             const register double tmp1_1 = X_0*w68;  
                             const register double tmp0_1 = X_1*w69;  
2106                              EM_F[0]+=tmp0_1 + tmp1_1;                              EM_F[0]+=tmp0_1 + tmp1_1;
2107                              EM_F[1]+=tmp0_1 + tmp2_1;                              EM_F[1]+=tmp0_1 + tmp2_1;
2108                              EM_F[2]+=tmp1_1 + tmp3_1;                              EM_F[2]+=tmp1_1 + tmp3_1;
# Line 2050  void Rectangle::assemblePDESingle(Paso_S Line 2122  void Rectangle::assemblePDESingle(Paso_S
2122                              const register double Y_3 = Y_p[3];                              const register double Y_3 = Y_p[3];
2123                              const register double tmp0_0 = Y_1 + Y_2;                              const register double tmp0_0 = Y_1 + Y_2;
2124                              const register double tmp1_0 = Y_0 + Y_3;                              const register double tmp1_0 = Y_0 + Y_3;
2125                              const register double tmp9_1 = Y_0*w74;                              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;                              const register double tmp4_1 = tmp1_0*w73;
2130                              const register double tmp5_1 = Y_2*w74;                              const register double tmp5_1 = Y_2*w74;
                             const register double tmp7_1 = Y_1*w74;  
2131                              const register double tmp6_1 = Y_2*w72;                              const register double tmp6_1 = Y_2*w72;
2132                              const register double tmp2_1 = Y_3*w74;                              const register double tmp7_1 = Y_1*w74;
2133                              const register double tmp8_1 = Y_3*w72;                              const register double tmp8_1 = Y_3*w72;
2134                              const register double tmp3_1 = Y_1*w72;                              const register double tmp9_1 = Y_0*w74;
                             const register double tmp0_1 = Y_0*w72;  
                             const register double tmp1_1 = tmp0_0*w73;  
2135                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2136                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2137                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2138                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2139                          } else { /* constant data */                          } else { // constant data
2140                              const register double Y_0 = Y_p[0];                              const register double tmp0_1 = Y_p[0]*w75;
                             const register double tmp0_1 = Y_0*w75;  
2141                              EM_F[0]+=tmp0_1;                              EM_F[0]+=tmp0_1;
2142                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
2143                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
2144                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
2145                          }                          }
2146                      }                      }
                     /* GENERATOR SNIP_PDE_SINGLE BOTTOM */  
2147    
2148                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2149                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
# Line 2083  void Rectangle::assemblePDESingle(Paso_S Line 2153  void Rectangle::assemblePDESingle(Paso_S
2153                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2154                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2155                      if (add_EM_F) {                      if (add_EM_F) {
                         //cout << "-- AddtoRHS -- " << endl;  
2156                          double *F_p=rhs.getSampleDataRW(0);                          double *F_p=rhs.getSampleDataRW(0);
2157                          for (index_t i=0; i<rowIndex.size(); i++) {                          for (index_t i=0; i<rowIndex.size(); i++) {
2158                              if (rowIndex[i]<getNumDOF()) {                              if (rowIndex[i]<getNumDOF()) {
2159                                  F_p[rowIndex[i]]+=EM_F[i];                                  F_p[rowIndex[i]]+=EM_F[i];
                                 //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;  
2160                              }                              }
2161                          }                          }
                         //cout << "---"<<endl;  
2162                      }                      }
2163                      if (add_EM_S) {                      if (add_EM_S) {
                         //cout << "-- AddtoSystem -- " << endl;  
2164                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2165                      }                      }
2166    
2167                  } // end k0 loop                  } // end k0 loop
2168              } // end k1 loop              } // end k1 loop
2169          } // end of coloring          } // end of colouring
2170      } // end of parallel region      } // 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;
2594                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
2595                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
2596                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
2597                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2598                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
2599                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
2600                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
2601                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
2602                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
2603                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
2604                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2605                                    }
2606                                }
2607                            } else { // constant data
2608                                for (index_t k=0; k<numEq; k++) {
2609                                    for (index_t m=0; m<numComp; m++) {
2610                                        const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
2611                                        const register double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
2612                                        const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
2613                                        const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
2614                                        const register double tmp0_0 = A_01 + A_10;
2615                                        const register double tmp0_1 = A_00*w18;
2616                                        const register double tmp1_1 = A_01*w19;
2617                                        const register double tmp2_1 = A_10*w20;
2618                                        const register double tmp3_1 = A_11*w21;
2619                                        const register double tmp4_1 = A_00*w22;
2620                                        const register double tmp5_1 = tmp0_0*w19;
2621                                        const register double tmp6_1 = A_11*w23;
2622                                        const register double tmp7_1 = A_11*w25;
2623                                        const register double tmp8_1 = A_00*w24;
2624                                        const register double tmp9_1 = tmp0_0*w20;
2625                                        const register double tmp10_1 = A_01*w20;
2626                                        const register double tmp11_1 = A_11*w27;
2627                                        const register double tmp12_1 = A_00*w26;
2628                                        const register double tmp13_1 = A_10*w19;
2629                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2630                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2631                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
2632                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2633                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2634                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2635                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
2636                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
2637                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
2638                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
2639                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2640                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2641                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2642                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
2643                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2644                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2645                                    }
2646                                }
2647                            }
2648                        }
2649                        ///////////////
2650                        // process B //
2651                        ///////////////
2652                        if (!B.isEmpty()) {
2653                            add_EM_S=true;
2654                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2655                            if (B.actsExpanded()) {
2656                                for (index_t k=0; k<numEq; k++) {
2657                                    for (index_t m=0; m<numComp; m++) {
2658                                        const register double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
2659                                        const register double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
2660                                        const register double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
2661                                        const register double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
2662                                        const register double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];
2663                                        const register double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
2664                                        const register double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
2665                                        const register double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
2666                                        const register double tmp0_0 = B_1_0 + B_1_1;
2667                                        const register double tmp1_0 = B_1_2 + B_1_3;
2668                                        const register double tmp2_0 = B_0_1 + B_0_3;
2669                                        const register double tmp3_0 = B_0_0 + B_0_2;
2670                                        const register double tmp63_1 = B_1_1*w42;
2671                                        const register double tmp79_1 = B_1_1*w40;
2672                                        const register double tmp37_1 = tmp3_0*w35;
2673                                        const register double tmp8_1 = tmp0_0*w32;
2674                                        const register double tmp71_1 = B_0_1*w34;
2675                                        const register double tmp19_1 = B_0_3*w31;
2676                                        const register double tmp15_1 = B_0_3*w34;
2677                                        const register double tmp9_1 = tmp3_0*w34;
2678                                        const register double tmp35_1 = B_1_0*w36;
2679                                        const register double tmp66_1 = B_0_3*w28;
2680                                        const register double tmp28_1 = B_1_0*w42;
2681                                        const register double tmp22_1 = B_1_0*w40;
2682                                        const register double tmp16_1 = B_1_2*w29;
2683                                        const register double tmp6_1 = tmp2_0*w35;
2684                                        const register double tmp55_1 = B_1_3*w40;
2685                                        const register double tmp50_1 = B_1_3*w42;
2686                                        const register double tmp7_1 = tmp1_0*w29;
2687                                        const register double tmp1_1 = tmp1_0*w32;
2688                                        const register double tmp57_1 = B_0_3*w30;
2689                                        const register double tmp18_1 = B_1_1*w32;
2690                                        const register double tmp53_1 = B_1_0*w41;
2691                                        const register double tmp61_1 = B_1_3*w36;
2692                                        const register double tmp27_1 = B_0_3*w38;
2693                                        const register double tmp64_1 = B_0_2*w30;
2694                                        const register double tmp76_1 = B_0_1*w38;
2695                                        const register double tmp39_1 = tmp2_0*w34;
2696                                        const register double tmp62_1 = B_0_1*w31;
2697                                        const register double tmp56_1 = B_0_0*w31;
2698                                        const register double tmp49_1 = B_1_1*w36;
2699                                        const register double tmp2_1 = B_0_2*w31;
2700                                        const register double tmp23_1 = B_0_2*w33;
2701                                        const register double tmp38_1 = B_1_1*w43;
2702                                        const register double tmp74_1 = B_1_2*w41;
2703                                        const register double tmp43_1 = B_1_1*w41;
2704                                        const register double tmp58_1 = B_0_2*w28;
2705                                        const register double tmp67_1 = B_0_0*w33;
2706                                        const register double tmp33_1 = tmp0_0*w39;
2707                                        const register double tmp4_1 = B_0_0*w28;
2708                                        const register double tmp20_1 = B_0_0*w30;
2709                                        const register double tmp13_1 = B_0_2*w38;
2710                                        const register double tmp65_1 = B_1_2*w43;
2711                                        const register double tmp0_1 = tmp0_0*w29;
2712                                        const register double tmp41_1 = tmp3_0*w33;
2713                                        const register double tmp73_1 = B_0_2*w37;
2714                                        const register double tmp69_1 = B_0_0*w38;
2715                                        const register double tmp48_1 = B_1_2*w39;
2716                                        const register double tmp59_1 = B_0_1*w33;
2717                                        const register double tmp17_1 = B_1_3*w41;
2718                                        const register double tmp5_1 = B_0_3*w33;
2719                                        const register double tmp3_1 = B_0_1*w30;
2720                                        const register double tmp21_1 = B_0_1*w28;
2721                                        const register double tmp42_1 = B_1_0*w29;
2722                                        const register double tmp54_1 = B_1_2*w32;
2723                                        const register double tmp60_1 = B_1_0*w39;
2724                                        const register double tmp32_1 = tmp1_0*w36;
2725                                        const register double tmp10_1 = B_0_1*w37;
2726                                        const register double tmp14_1 = B_0_0*w35;
2727                                        const register double tmp29_1 = B_0_1*w35;
2728                                        const register double tmp26_1 = B_1_2*w36;
2729                                        const register double tmp30_1 = B_1_3*w43;
2730                                        const register double tmp70_1 = B_0_2*w35;
2731                                        const register double tmp34_1 = B_1_3*w39;
2732                                        const register double tmp51_1 = B_1_0*w43;
2733                                        const register double tmp31_1 = B_0_2*w34;
2734                                        const register double tmp45_1 = tmp3_0*w28;
2735                                        const register double tmp11_1 = tmp1_0*w39;
2736                                        const register double tmp52_1 = B_1_1*w29;
2737                                        const register double tmp44_1 = B_1_3*w32;
2738                                        const register double tmp25_1 = B_1_1*w39;
2739                                        const register double tmp47_1 = tmp2_0*w33;
2740                                        const register double tmp72_1 = B_1_3*w29;
2741                                        const register double tmp40_1 = tmp2_0*w28;
2742                                        const register double tmp46_1 = B_1_2*w40;
2743                                        const register double tmp36_1 = B_1_2*w42;
2744                                        const register double tmp24_1 = B_0_0*w37;
2745                                        const register double tmp77_1 = B_0_3*w35;
2746                                        const register double tmp68_1 = B_0_3*w37;
2747                                        const register double tmp78_1 = B_0_0*w34;
2748                                        const register double tmp12_1 = tmp0_0*w36;
2749                                        const register double tmp75_1 = B_1_0*w32;
2750                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2751                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2752                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2753                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
2754                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2755                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
2756                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
2757                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
2758                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
2759                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2760                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
2761                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2762                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
2763                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2764                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
2765                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
2766                                    }
2767                                }
2768                            } else { // constant data
2769                                for (index_t k=0; k<numEq; k++) {
2770                                    for (index_t m=0; m<numComp; m++) {
2771                                        const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
2772                                        const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
2773                                        const register double tmp6_1 = B_1*w50;
2774                                        const register double tmp1_1 = B_1*w45;
2775                                        const register double tmp5_1 = B_1*w49;
2776                                        const register double tmp4_1 = B_1*w48;
2777                                        const register double tmp0_1 = B_0*w44;
2778                                        const register double tmp2_1 = B_0*w46;
2779                                        const register double tmp7_1 = B_0*w51;
2780                                        const register double tmp3_1 = B_0*w47;
2781                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
2782                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
2783                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
2784                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;
2785                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;
2786                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
2787                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;
2788                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;
2789                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;
2790                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;
2791                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;
2792                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;
2793                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;
2794                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
2795                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;
2796                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;
2797                                    }
2798                                }
2799                            }
2800                        }
2801                        ///////////////
2802                        // process C //
2803                        ///////////////
2804                        if (!C.isEmpty()) {
2805                            add_EM_S=true;
2806                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2807                            if (C.actsExpanded()) {
2808                                for (index_t k=0; k<numEq; k++) {
2809                                    for (index_t m=0; m<numComp; m++) {
2810                                        const register double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
2811                                        const register double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
2812                                        const register double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
2813                                        const register double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
2814                                        const register double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];
2815                                        const register double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
2816                                        const register double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];
2817                                        const register double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];
2818                                        const register double tmp2_0 = C_0_1 + C_0_3;
2819                                        const register double tmp1_0 = C_1_2 + C_1_3;
2820                                        const register double tmp3_0 = C_0_0 + C_0_2;
2821                                        const register double tmp0_0 = C_1_0 + C_1_1;
2822                                        const register double tmp64_1 = C_0_2*w30;
2823                                        const register double tmp14_1 = C_0_2*w28;
2824                                        const register double tmp19_1 = C_0_3*w31;
2825                                        const register double tmp22_1 = C_1_0*w40;
2826                                        const register double tmp37_1 = tmp3_0*w35;
2827                                        const register double tmp29_1 = C_0_1*w35;
2828                                        const register double tmp73_1 = C_0_2*w37;
2829                                        const register double tmp74_1 = C_1_2*w41;
2830                                        const register double tmp52_1 = C_1_3*w39;
2831                                        const register double tmp25_1 = C_1_1*w39;
2832                                        const register double tmp62_1 = C_0_1*w31;
2833                                        const register double tmp79_1 = C_1_1*w40;
2834                                        const register double tmp43_1 = C_1_1*w36;
2835                                        const register double tmp27_1 = C_0_3*w38;
2836                                        const register double tmp28_1 = C_1_0*w42;
2837                                        const register double tmp63_1 = C_1_1*w42;
2838                                        const register double tmp59_1 = C_0_3*w34;
2839                                        const register double tmp72_1 = C_1_3*w29;
2840                                        const register double tmp40_1 = tmp2_0*w35;
2841                                        const register double tmp13_1 = C_0_3*w30;
2842                                        const register double tmp51_1 = C_1_2*w40;
2843                                        const register double tmp54_1 = C_1_2*w42;
2844                                        const register double tmp12_1 = C_0_0*w31;
2845                                        const register double tmp2_1 = tmp1_0*w32;
2846                                        const register double tmp68_1 = C_0_2*w31;
2847                                        const register double tmp75_1 = C_1_0*w32;
2848                                        const register double tmp49_1 = C_1_1*w41;
2849                                        const register double tmp4_1 = C_0_2*w35;
2850                                        const register double tmp66_1 = C_0_3*w28;
2851                                        const register double tmp56_1 = C_0_1*w37;
2852                                        const register double tmp5_1 = C_0_1*w34;
2853                                        const register double tmp38_1 = tmp2_0*w34;
2854                                        const register double tmp76_1 = C_0_1*w38;
2855                                        const register double tmp21_1 = C_0_1*w28;
2856                                        const register double tmp69_1 = C_0_1*w30;
2857                                        const register double tmp53_1 = C_1_0*w36;
2858                                        const register double tmp42_1 = C_1_2*w39;
2859                                        const register double tmp32_1 = tmp1_0*w29;
2860                                        const register double tmp45_1 = C_1_0*w43;
2861                                        const register double tmp33_1 = tmp0_0*w32;
2862                                        const register double tmp35_1 = C_1_0*w41;
2863                                        const register double tmp26_1 = C_1_2*w36;
2864                                        const register double tmp67_1 = C_0_0*w33;
2865                                        const register double tmp31_1 = C_0_2*w34;
2866                                        const register double tmp20_1 = C_0_0*w30;
2867                                        const register double tmp70_1 = C_0_0*w28;
2868                                        const register double tmp8_1 = tmp0_0*w39;
2869                                        const register double tmp30_1 = C_1_3*w43;
2870                                        const register double tmp0_1 = tmp0_0*w29;
2871                                        const register double tmp17_1 = C_1_3*w41;
2872                                        const register double tmp58_1 = C_0_0*w35;
2873                                        const register double tmp9_1 = tmp3_0*w33;
2874                                        const register double tmp61_1 = C_1_3*w36;
2875                                        const register double tmp41_1 = tmp3_0*w34;
2876                                        const register double tmp50_1 = C_1_3*w32;
2877                                        const register double tmp18_1 = C_1_1*w32;
2878                                        const register double tmp6_1 = tmp1_0*w36;
2879                                        const register double tmp3_1 = C_0_0*w38;
2880                                        const register double tmp34_1 = C_1_1*w29;
2881                                        const register double tmp77_1 = C_0_3*w35;
2882                                        const register double tmp65_1 = C_1_2*w43;
2883                                        const register double tmp71_1 = C_0_3*w33;
2884                                        const register double tmp55_1 = C_1_1*w43;
2885                                        const register double tmp46_1 = tmp3_0*w28;
2886                                        const register double tmp24_1 = C_0_0*w37;
2887                                        const register double tmp10_1 = tmp1_0*w39;
2888                                        const register double tmp48_1 = C_1_0*w29;
2889                                        const register double tmp15_1 = C_0_1*w33;
2890                                        const register double tmp36_1 = C_1_2*w32;
2891                                        const register double tmp60_1 = C_1_0*w39;
2892                                        const register double tmp47_1 = tmp2_0*w33;
2893                                        const register double tmp16_1 = C_1_2*w29;
2894                                        const register double tmp1_1 = C_0_3*w37;
2895                                        const register double tmp7_1 = tmp2_0*w28;
2896                                        const register double tmp39_1 = C_1_3*w40;
2897                                        const register double tmp44_1 = C_1_3*w42;
2898                                        const register double tmp57_1 = C_0_2*w38;
2899                                        const register double tmp78_1 = C_0_0*w34;
2900                                        const register double tmp11_1 = tmp0_0*w36;
2901                                        const register double tmp23_1 = C_0_2*w33;
2902                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2903                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2904                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2905                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
2906                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2907                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
2908                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
2909                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
2910                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
2911                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2912                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
2913                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2914                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
2915                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2916                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
2917                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
2918                                    }
2919                                }
2920                            } else { // constant data
2921                                for (index_t k=0; k<numEq; k++) {
2922                                    for (index_t m=0; m<numComp; m++) {
2923                                        const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
2924                                        const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
2925                                        const register double tmp1_1 = C_1*w45;
2926                                        const register double tmp3_1 = C_0*w51;
2927                                        const register double tmp4_1 = C_0*w44;
2928                                        const register double tmp7_1 = C_0*w46;
2929                                        const register double tmp5_1 = C_1*w49;
2930                                        const register double tmp2_1 = C_1*w48;
2931                                        const register double tmp0_1 = C_0*w47;
2932                                        const register double tmp6_1 = C_1*w50;
2933                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
2934                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
2935                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
2936                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp4_1 + tmp5_1;
2937                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp6_1;
2938                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
2939                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp7_1;
2940                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp7_1;
2941                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp3_1 + tmp6_1;
2942                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp3_1 + tmp5_1;
2943                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp6_1 + tmp7_1;
2944                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
2945                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp4_1 + tmp6_1;
2946                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp4_1;
2947                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2_1 + tmp7_1;
2948                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp5_1;
2949                                    }
2950                                }
2951                            }
2952                        }
2953                        ///////////////
2954                        // process D //
2955                        ///////////////
2956                        if (!D.isEmpty()) {
2957                            add_EM_S=true;
2958                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2959                            if (D.actsExpanded()) {
2960                                for (index_t k=0; k<numEq; k++) {
2961                                    for (index_t m=0; m<numComp; m++) {
2962                                        const register double D_0 = D_p[INDEX3(k, m, 0, numEq, numComp)];
2963                                        const register double D_1 = D_p[INDEX3(k, m, 1, numEq, numComp)];
2964                                        const register double D_2 = D_p[INDEX3(k, m, 2, numEq, numComp)];
2965                                        const register double D_3 = D_p[INDEX3(k, m, 3, numEq, numComp)];
2966                                        const register double tmp4_0 = D_1 + D_3;
2967                                        const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;
2968                                        const register double tmp5_0 = D_0 + D_2;
2969                                        const register double tmp0_0 = D_0 + D_1;
2970                                        const register double tmp6_0 = D_0 + D_3;
2971                                        const register double tmp1_0 = D_2 + D_3;
2972                                        const register double tmp3_0 = D_1 + D_2;
2973                                        const register double tmp16_1 = D_1*w56;
2974                                        const register double tmp14_1 = tmp6_0*w54;
2975                                        const register double tmp8_1 = D_3*w55;
2976                                        const register double tmp2_1 = tmp2_0*w54;
2977                                        const register double tmp12_1 = tmp5_0*w52;
2978                                        const register double tmp4_1 = tmp0_0*w53;
2979                                        const register double tmp3_1 = tmp1_0*w52;
2980                                        const register double tmp13_1 = tmp4_0*w53;
2981                                        const register double tmp10_1 = tmp4_0*w52;
2982                                        const register double tmp1_1 = tmp1_0*w53;
2983                                        const register double tmp7_1 = D_3*w56;
2984                                        const register double tmp0_1 = tmp0_0*w52;
2985                                        const register double tmp11_1 = tmp5_0*w53;
2986                                        const register double tmp9_1 = D_0*w56;
2987                                        const register double tmp5_1 = tmp3_0*w54;
2988                                        const register double tmp18_1 = D_2*w56;
2989                                        const register double tmp17_1 = D_1*w55;
2990                                        const register double tmp6_1 = D_0*w55;
2991                                        const register double tmp15_1 = D_2*w55;
2992                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
2993                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1;
2994                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
2995                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2996                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2997                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1;
2998                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1;
2999                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp2_1;
3000                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp12_1 + tmp13_1;
3001                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp12_1 + tmp13_1;
3002                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp10_1 + tmp11_1;
3003                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
3004                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
3005                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3006                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2_1;
3007                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
3008                                    }
3009                                 }
3010                            } else { // constant data
3011                                for (index_t k=0; k<numEq; k++) {
3012                                    for (index_t m=0; m<numComp; m++) {
3013                                        const register double D_0 = D_p[INDEX2(k, m, numEq)];
3014                                        const register double tmp2_1 = D_0*w59;
3015                                        const register double tmp1_1 = D_0*w58;
3016                                        const register double tmp0_1 = D_0*w57;
3017                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
3018                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1;
3019                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
3020                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp2_1;
3021                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1;
3022                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1;
3023                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;
3024                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1;
3025                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;
3026                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;
3027                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;
3028                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
3029                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1;
3030                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
3031                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1;
3032                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp2_1;
3033                                    }
3034                                }
3035                            }
3036                        }
3037                        ///////////////
3038                        // process X //
3039                        ///////////////
3040                        if (!X.isEmpty()) {
3041                            add_EM_F=true;
3042                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3043                            if (X.actsExpanded()) {
3044                                for (index_t k=0; k<numEq; k++) {
3045                                    const register double X_0_0 = X_p[INDEX3(k, 0, 0, numEq, 2)];
3046                                    const register double X_1_0 = X_p[INDEX3(k, 1, 0, numEq, 2)];
3047                                    const register double X_0_1 = X_p[INDEX3(k, 0, 1, numEq, 2)];
3048                                    const register double X_1_1 = X_p[INDEX3(k, 1, 1, numEq, 2)];
3049                                    const register double X_0_2 = X_p[INDEX3(k, 0, 2, numEq, 2)];
3050                                    const register double X_1_2 = X_p[INDEX3(k, 1, 2, numEq, 2)];
3051                                    const register double X_0_3 = X_p[INDEX3(k, 0, 3, numEq, 2)];
3052                                    const register double X_1_3 = X_p[INDEX3(k, 1, 3, numEq, 2)];
3053                                    const register double tmp1_0 = X_1_1 + X_1_3;
3054                                    const register double tmp3_0 = X_0_0 + X_0_1;
3055                                    const register double tmp2_0 = X_1_0 + X_1_2;
3056                                    const register double tmp0_0 = X_0_2 + X_0_3;
3057                                    const register double tmp8_1 = tmp2_0*w66;
3058                                    const register double tmp5_1 = tmp3_0*w64;
3059                                    const register double tmp14_1 = tmp0_0*w64;
3060                                    const register double tmp3_1 = tmp3_0*w60;
3061                                    const register double tmp9_1 = tmp3_0*w63;
3062                                    const register double tmp13_1 = tmp3_0*w65;
3063                                    const register double tmp12_1 = tmp1_0*w66;
3064                                    const register double tmp10_1 = tmp0_0*w60;
3065                                    const register double tmp2_1 = tmp2_0*w61;
3066                                    const register double tmp6_1 = tmp2_0*w62;
3067                                    const register double tmp4_1 = tmp0_0*w65;
3068                                    const register double tmp11_1 = tmp1_0*w67;
3069                                    const register double tmp1_1 = tmp1_0*w62;
3070                                    const register double tmp7_1 = tmp1_0*w61;
3071                                    const register double tmp0_1 = tmp0_0*w63;
3072                                    const register double tmp15_1 = tmp2_0*w67;
3073                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3074                                    EM_F[INDEX2(k,1,numEq)]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
3075                                    EM_F[INDEX2(k,2,numEq)]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
3076                                    EM_F[INDEX2(k,3,numEq)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
3077                                }
3078                            } else { // constant data
3079                                for (index_t k=0; k<numEq; k++) {
3080                                    const register double X_0 = X_p[INDEX2(k, 0, numEq)];
3081                                    const register double X_1 = X_p[INDEX2(k, 1, numEq)];
3082                                    const register double tmp0_1 = X_1*w69;
3083                                    const register double tmp1_1 = X_0*w68;
3084                                    const register double tmp2_1 = X_0*w70;
3085                                    const register double tmp3_1 = X_1*w71;
3086                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3087                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;
3088                                    EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;
3089                                    EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;
3090                                }
3091                            }
3092                        }
3093                        ///////////////
3094                        // process Y //
3095                        ///////////////
3096                        if (!Y.isEmpty()) {
3097                            add_EM_F=true;
3098                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3099                            if (Y.actsExpanded()) {
3100                                for (index_t k=0; k<numEq; k++) {
3101                                    const register double Y_0 = Y_p[INDEX2(k, 0, numEq)];
3102                                    const register double Y_1 = Y_p[INDEX2(k, 1, numEq)];
3103                                    const register double Y_2 = Y_p[INDEX2(k, 2, numEq)];
3104                                    const register double Y_3 = Y_p[INDEX2(k, 3, numEq)];
3105                                    const register double tmp0_0 = Y_1 + Y_2;
3106                                    const register double tmp1_0 = Y_0 + Y_3;
3107                                    const register double tmp0_1 = Y_0*w72;
3108                                    const register double tmp1_1 = tmp0_0*w73;
3109                                    const register double tmp2_1 = Y_3*w74;
3110                                    const register double tmp3_1 = Y_1*w72;
3111                                    const register double tmp4_1 = tmp1_0*w73;
3112                                    const register double tmp5_1 = Y_2*w74;
3113                                    const register double tmp6_1 = Y_2*w72;
3114                                    const register double tmp7_1 = Y_1*w74;
3115                                    const register double tmp8_1 = Y_3*w72;
3116                                    const register double tmp9_1 = Y_0*w74;
3117                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
3118                                    EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
3119                                    EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
3120                                    EM_F[INDEX2(k,3,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;
3121                                }
3122                            } else { // constant data
3123                                for (index_t k=0; k<numEq; k++) {
3124                                    const register double tmp0_1 = Y_p[k]*w75;
3125                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3126                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3127                                    EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
3128                                    EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
3129                                }
3130                            }
3131                        }
3132    
3133                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
3134                        const index_t firstNode=m_N0*k1+k0;
3135                        IndexVector rowIndex;
3136                        rowIndex.push_back(m_dofMap[firstNode]);
3137                        rowIndex.push_back(m_dofMap[firstNode+1]);
3138                        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
3139                        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
3140                        if (add_EM_F) {
3141                            double *F_p=rhs.getSampleDataRW(0);
3142                            for (index_t i=0; i<rowIndex.size(); i++) {
3143                                if (rowIndex[i]<getNumDOF()) {
3144                                    for (index_t eq=0; eq<numEq; eq++) {
3145                                        F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];
3146                                    }
3147                                }
3148                            }
3149                        }
3150                        if (add_EM_S) {
3151                            addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);
3152                        }
3153    
3154                    } // end k0 loop
3155                } // end k1 loop
3156            } // end of colouring
3157        } // end of parallel region
3158    }
3159    
3160    //protected
3161    void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,
3162            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
3163            const escript::Data& C, const escript::Data& D,
3164            const escript::Data& X, const escript::Data& Y,
3165            const escript::Data& d, const escript::Data& y) const
3166    {
3167        const double h0 = m_l0/m_gNE0;
3168        const double h1 = m_l1/m_gNE1;
3169        dim_t numEq, numComp;
3170        if (!mat)
3171            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
3172        else {
3173            numEq=mat->logical_row_block_size;
3174            numComp=mat->logical_col_block_size;
3175        }
3176    
3177        const double w0 = -.25*h1/h0;
3178        const double w1 = .25;
3179        const double w2 = -.25;
3180        const double w3 = .25*h0/h1;
3181        const double w4 = -.25*h0/h1;
3182        const double w5 = .25*h1/h0;
3183        const double w6 = -.125*h1;
3184        const double w7 = -.125*h0;
3185        const double w8 = .125*h1;
3186        const double w9 = .125*h0;
3187        const double w10 = .0625*h0*h1;
3188        const double w11 = -.5*h1;
3189        const double w12 = -.5*h0;
3190        const double w13 = .5*h1;
3191        const double w14 = .5*h0;
3192        const double w15 = .25*h0*h1;
3193    
3194        rhs.requireWrite();
3195    #pragma omp parallel
3196        {
3197            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3198    #pragma omp for
3199                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3200                    for (index_t k0=0; k0<m_NE0; ++k0)  {
3201                        bool add_EM_S=false;
3202                        bool add_EM_F=false;
3203                        vector<double> EM_S(4*4*numEq*numComp, 0);
3204                        vector<double> EM_F(4*numEq, 0);
3205                        const index_t e = k0 + m_NE0*k1;
3206                        ///////////////
3207                        // process A //
3208                        ///////////////
3209                        if (!A.isEmpty()) {
3210                            add_EM_S=true;
3211                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
3212                            for (index_t k=0; k<numEq; k++) {
3213                                for (index_t m=0; m<numComp; m++) {
3214                                    const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
3215                                    const register double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
3216                                    const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
3217                                    const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
3218                                    const register double tmp0_0 = A_01 + A_10;
3219                                    const register double tmp0_1 = A_11*w3;
3220                                    const register double tmp1_1 = A_00*w0;
3221                                    const register double tmp2_1 = A_01*w1;
3222                                    const register double tmp3_1 = A_10*w2;
3223                                    const register double tmp4_1 = tmp0_0*w1;
3224                                    const register double tmp5_1 = A_11*w4;
3225                                    const register double tmp6_1 = A_00*w5;
3226                                    const register double tmp7_1 = tmp0_0*w2;
3227                                    const register double tmp8_1 = A_10*w1;
3228                                    const register double tmp9_1 = A_01*w2;
3229                                    EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
3230                                    EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
3231                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
3232                                    EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
3233                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3234                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
3235                                    EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
3236                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
3237                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
3238                                    EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
3239                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
3240                                    EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3241                                    EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
3242                                    EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
3243                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
3244                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
3245                                }
3246                            }
3247                        }
3248                        ///////////////
3249                        // process B //
3250                        ///////////////
3251                        if (!B.isEmpty()) {
3252                            add_EM_S=true;
3253                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
3254                            for (index_t k=0; k<numEq; k++) {
3255                                for (index_t m=0; m<numComp; m++) {
3256                                    const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
3257                                    const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
3258                                    const register double tmp0_1 = B_0*w6;
3259                                    const register double tmp1_1 = B_1*w7;
3260                                    const register double tmp2_1 = B_0*w8;
3261                                    const register double tmp3_1 = B_1*w9;
3262                                    EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3263                                    EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3264                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3265                                    EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3266                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3267                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3268                                    EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3269                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3270                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3271                                    EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3272                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3273                                    EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3274                                    EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3275                                    EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3276                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3277                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3278                                }
3279                            }
3280                        }
3281                        ///////////////
3282                        // process C //
3283                        ///////////////
3284                        if (!C.isEmpty()) {
3285                            add_EM_S=true;
3286                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
3287                            for (index_t k=0; k<numEq; k++) {
3288                                for (index_t m=0; m<numComp; m++) {
3289                                    const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
3290                                    const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
3291                                    const register double tmp0_1 = C_0*w8;
3292                                    const register double tmp1_1 = C_1*w7;
3293                                    const register double tmp2_1 = C_1*w9;
3294                                    const register double tmp3_1 = C_0*w6;
3295                                    EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3296                                    EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3297                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3298                                    EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3299                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3300                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3301                                    EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3302                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3303                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3304                                    EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3305                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3306                                    EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3307                                    EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3308                                    EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3309                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3310                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3311                                }
3312                            }
3313                        }
3314                        ///////////////
3315                        // process D //
3316                        ///////////////
3317                        if (!D.isEmpty()) {
3318                            add_EM_S=true;
3319                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
3320                            for (index_t k=0; k<numEq; k++) {
3321                                for (index_t m=0; m<numComp; m++) {
3322                                    const register double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;
3323                                    EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
3324                                    EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
3325                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;
3326                                    EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;
3327                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
3328                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
3329                                    EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1;
3330                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;
3331                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;
3332                                    EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;
3333                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
3334                                    EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
3335                                    EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1;
3336                                    EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;
3337                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
3338                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
3339                                }
3340                            }
3341                        }
3342                        ///////////////
3343                        // process X //
3344                        ///////////////
3345                        if (!X.isEmpty()) {
3346                            add_EM_F=true;
3347                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3348                            for (index_t k=0; k<numEq; k++) {
3349                                const register double X_0 = X_p[INDEX2(k, 0, numEq)];
3350                                const register double X_1 = X_p[INDEX2(k, 1, numEq)];
3351                                const register double tmp0_1 = X_0*w11;
3352                                const register double tmp1_1 = X_1*w12;
3353                                const register double tmp2_1 = X_0*w13;
3354                                const register double tmp3_1 = X_1*w14;
3355                                EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3356                                EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;
3357                                EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp3_1;
3358                                EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;
3359                            }
3360                        }
3361                        ///////////////
3362                        // process Y //
3363                        ///////////////
3364                        if (!Y.isEmpty()) {
3365                            add_EM_F=true;
3366                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3367                            for (index_t k=0; k<numEq; k++) {
3368                                const register double tmp0_1 = Y_p[k]*w15;
3369                                EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3370                                EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3371                                EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
3372         &nb