/[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 3756 by caltinay, Fri Jan 6 02:35:19 2012 UTC revision 3777 by caltinay, Thu Jan 19 06:17:38 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 106  void Rectangle::dump(const string& fileN Line 111  void Rectangle::dump(const string& fileN
111      const int NUM_SILO_FILES = 1;      const int NUM_SILO_FILES = 1;
112      const char* blockDirFmt = "/block%04d";      const char* blockDirFmt = "/block%04d";
113      int driver=DB_HDF5;          int driver=DB_HDF5;    
     string siloPath;  
114      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
115    
116  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 126  void Rectangle::dump(const string& fileN Line 130  void Rectangle::dump(const string& fileN
130                          PMPIO_DefaultClose, (void*)&driver);                          PMPIO_DefaultClose, (void*)&driver);
131          }          }
132          if (baton) {          if (baton) {
133              char str[64];              char siloPath[64];
134              snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));              snprintf(siloPath, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
135              siloPath = str;              dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath);
             dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());  
136          }          }
137  #endif  #endif
138      } else {      } else {
# Line 141  void Rectangle::dump(const string& fileN Line 144  void Rectangle::dump(const string& fileN
144              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
145                      getDescription().c_str(), driver);                      getDescription().c_str(), driver);
146          }          }
147            char siloPath[64];
148            snprintf(siloPath, 64, blockDirFmt, 0);
149            DBMkDir(dbfile, siloPath);
150            DBSetDir(dbfile, siloPath);
151      }      }
152    
153      if (!dbfile)      if (!dbfile)
# Line 241  const int* Rectangle::borrowSampleRefere Line 248  const int* Rectangle::borrowSampleRefere
248  {  {
249      switch (fsType) {      switch (fsType) {
250          case Nodes:          case Nodes:
251          case ReducedNodes: //FIXME: reduced          case ReducedNodes: // FIXME: reduced
252              return &m_nodeId[0];              return &m_nodeId[0];
253          case DegreesOfFreedom:          case DegreesOfFreedom:
254          case ReducedDegreesOfFreedom: //FIXME: reduced          case ReducedDegreesOfFreedom: // FIXME: reduced
255              return &m_dofId[0];              return &m_dofId[0];
256          case Elements:          case Elements:
257          case ReducedElements:          case ReducedElements:
# Line 262  const int* Rectangle::borrowSampleRefere Line 269  const int* Rectangle::borrowSampleRefere
269      throw RipleyException(msg.str());      throw RipleyException(msg.str());
270  }  }
271    
272  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsType, index_t id) const
273  {  {
274  #ifdef ESYS_MPI      if (getMPISize()==1)
275      if (fsCode == Nodes) {          return true;
276          return (m_dofMap[id] < getNumDOF());  
277        switch (fsType) {
278            case Nodes:
279            case ReducedNodes: // FIXME: reduced
280                return (m_dofMap[id] < getNumDOF());
281            case DegreesOfFreedom:
282            case ReducedDegreesOfFreedom:
283                return true;
284            case Elements:
285            case ReducedElements:
286                // check ownership of element's bottom left node
287                return (m_dofMap[id%m_NE0+m_N0*(id/m_NE0)] < getNumDOF());
288            case FaceElements:
289            case ReducedFaceElements:
290                {
291                    // determine which face the sample belongs to before
292                    // checking ownership of corresponding element's first node
293                    const IndexVector faces = getNumFacesPerBoundary();
294                    dim_t n=0;
295                    for (size_t i=0; i<faces.size(); i++) {
296                        n+=faces[i];
297                        if (id<n) {
298                            index_t k;
299                            if (i==1)
300                                k=m_N0-2;
301                            else if (i==3)
302                                k=m_N0*(m_N1-2);
303                            else
304                                k=0;
305                            // determine whether to move right or up
306                            const index_t delta=(i/2==0 ? m_N0 : 1);
307                            return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());
308                        }
309                    }
310                    return false;
311                }
312            default:
313                break;
314        }
315    
316        stringstream msg;
317        msg << "ownSample() not implemented for "
318            << functionSpaceTypeAsString(fsType);
319        throw RipleyException(msg.str());
320    }
321    
322    void Rectangle::setToNormal(escript::Data& out) const
323    {
324        if (out.getFunctionSpace().getTypeCode() == FaceElements) {
325            out.requireWrite();
326    #pragma omp parallel
327            {
328                if (m_faceOffset[0] > -1) {
329    #pragma omp for nowait
330                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
331                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
332                        // set vector at two quadrature points
333                        *o++ = -1.;
334                        *o++ = 0.;
335                        *o++ = -1.;
336                        *o = 0.;
337                    }
338                }
339    
340                if (m_faceOffset[1] > -1) {
341    #pragma omp for nowait
342                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
343                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
344                        // set vector at two quadrature points
345                        *o++ = 1.;
346                        *o++ = 0.;
347                        *o++ = 1.;
348                        *o = 0.;
349                    }
350                }
351    
352                if (m_faceOffset[2] > -1) {
353    #pragma omp for nowait
354                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
355                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
356                        // set vector at two quadrature points
357                        *o++ = 0.;
358                        *o++ = -1.;
359                        *o++ = 0.;
360                        *o = -1.;
361                    }
362                }
363    
364                if (m_faceOffset[3] > -1) {
365    #pragma omp for nowait
366                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
367                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
368                        // set vector at two quadrature points
369                        *o++ = 0.;
370                        *o++ = 1.;
371                        *o++ = 0.;
372                        *o = 1.;
373                    }
374                }
375            } // end of parallel section
376        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
377            out.requireWrite();
378    #pragma omp parallel
379            {
380                if (m_faceOffset[0] > -1) {
381    #pragma omp for nowait
382                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
383                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
384                        *o++ = -1.;
385                        *o = 0.;
386                    }
387                }
388    
389                if (m_faceOffset[1] > -1) {
390    #pragma omp for nowait
391                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
392                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
393                        *o++ = 1.;
394                        *o = 0.;
395                    }
396                }
397    
398                if (m_faceOffset[2] > -1) {
399    #pragma omp for nowait
400                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
401                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
402                        *o++ = 0.;
403                        *o = -1.;
404                    }
405                }
406    
407                if (m_faceOffset[3] > -1) {
408    #pragma omp for nowait
409                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
410                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
411                        *o++ = 0.;
412                        *o = 1.;
413                    }
414                }
415            } // end of parallel section
416    
417      } else {      } else {
418          stringstream msg;          stringstream msg;
419          msg << "ownSample() not implemented for "          msg << "setToNormal() not implemented for "
420              << functionSpaceTypeAsString(fsCode);              << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
421          throw RipleyException(msg.str());          throw RipleyException(msg.str());
422      }      }
 #else  
     return true;  
 #endif  
423  }  }
424    
425  void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const  void Rectangle::setToSize(escript::Data& out) const
426    {
427        if (out.getFunctionSpace().getTypeCode() == Elements
428                || out.getFunctionSpace().getTypeCode() == ReducedElements) {
429            out.requireWrite();
430            const dim_t numQuad=out.getNumDataPointsPerSample();
431            const double hSize=getFirstCoordAndSpacing(0).second;
432            const double vSize=getFirstCoordAndSpacing(1).second;
433            const double size=min(hSize,vSize);
434    #pragma omp parallel for
435            for (index_t k = 0; k < getNumElements(); ++k) {
436                double* o = out.getSampleDataRW(k);
437                fill(o, o+numQuad, size);
438            }
439        } else if (out.getFunctionSpace().getTypeCode() == FaceElements
440                || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
441            out.requireWrite();
442            const dim_t numQuad=out.getNumDataPointsPerSample();
443            const double hSize=getFirstCoordAndSpacing(0).second;
444            const double vSize=getFirstCoordAndSpacing(1).second;
445    #pragma omp parallel
446            {
447                if (m_faceOffset[0] > -1) {
448    #pragma omp for nowait
449                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
450                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
451                        fill(o, o+numQuad, vSize);
452                    }
453                }
454    
455                if (m_faceOffset[1] > -1) {
456    #pragma omp for nowait
457                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
458                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
459                        fill(o, o+numQuad, vSize);
460                    }
461                }
462    
463                if (m_faceOffset[2] > -1) {
464    #pragma omp for nowait
465                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
466                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
467                        fill(o, o+numQuad, hSize);
468                    }
469                }
470    
471                if (m_faceOffset[3] > -1) {
472    #pragma omp for nowait
473                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
474                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
475                        fill(o, o+numQuad, hSize);
476                    }
477                }
478            } // end of parallel section
479    
480        } else {
481            stringstream msg;
482            msg << "setToSize() not implemented for "
483                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
484            throw RipleyException(msg.str());
485        }
486    }
487    
488    Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
489                                                    bool reducedColOrder) const
490    {
491        /* FIXME: reduced
492        if (reducedRowOrder || reducedColOrder)
493            throw RipleyException("getPattern() not implemented for reduced order");
494        */
495        return m_pattern;
496    }
497    
498    void Rectangle::Print_Mesh_Info(const bool full) const
499    {
500        RipleyDomain::Print_Mesh_Info(full);
501        if (full) {
502            cout << "     Id  Coordinates" << endl;
503            cout.precision(15);
504            cout.setf(ios::scientific, ios::floatfield);
505            pair<double,double> xdx = getFirstCoordAndSpacing(0);
506            pair<double,double> ydy = getFirstCoordAndSpacing(1);
507            for (index_t i=0; i < getNumNodes(); i++) {
508                cout << "  " << setw(5) << m_nodeId[i]
509                    << "  " << xdx.first+(i%m_N0)*xdx.second
510                    << "  " << ydy.first+(i/m_N0)*ydy.second << endl;
511            }
512        }
513    }
514    
515    IndexVector Rectangle::getNumNodesPerDim() const
516    {
517        IndexVector ret;
518        ret.push_back(m_N0);
519        ret.push_back(m_N1);
520        return ret;
521    }
522    
523    IndexVector Rectangle::getNumElementsPerDim() const
524    {
525        IndexVector ret;
526        ret.push_back(m_NE0);
527        ret.push_back(m_NE1);
528        return ret;
529    }
530    
531    IndexVector Rectangle::getNumFacesPerBoundary() const
532    {
533        IndexVector ret(4, 0);
534        //left
535        if (m_offset0==0)
536            ret[0]=m_NE1;
537        //right
538        if (m_mpiInfo->rank%m_NX==m_NX-1)
539            ret[1]=m_NE1;
540        //bottom
541        if (m_offset1==0)
542            ret[2]=m_NE0;
543        //top
544        if (m_mpiInfo->rank/m_NX==m_NY-1)
545            ret[3]=m_NE0;
546        return ret;
547    }
548    
549    IndexVector Rectangle::getNumSubdivisionsPerDim() const
550    {
551        IndexVector ret;
552        ret.push_back(m_NX);
553        ret.push_back(m_NY);
554        return ret;
555    }
556    
557    pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
558    {
559        if (dim==0) {
560            return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
561        } else if (dim==1) {
562            return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
563        }
564        throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
565    }
566    
567    //protected
568    dim_t Rectangle::getNumDOF() const
569    {
570        return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;
571    }
572    
573    //protected
574    dim_t Rectangle::getNumFaceElements() const
575    {
576        const IndexVector faces = getNumFacesPerBoundary();
577        dim_t n=0;
578        for (size_t i=0; i<faces.size(); i++)
579            n+=faces[i];
580        return n;
581    }
582    
583    //protected
584    void Rectangle::assembleCoordinates(escript::Data& arg) const
585    {
586        escriptDataC x = arg.getDataC();
587        int numDim = m_numDim;
588        if (!isDataPointShapeEqual(&x, 1, &numDim))
589            throw RipleyException("setToX: Invalid Data object shape");
590        if (!numSamplesEqual(&x, 1, getNumNodes()))
591            throw RipleyException("setToX: Illegal number of samples in Data object");
592    
593        pair<double,double> xdx = getFirstCoordAndSpacing(0);
594        pair<double,double> ydy = getFirstCoordAndSpacing(1);
595        arg.requireWrite();
596    #pragma omp parallel for
597        for (dim_t i1 = 0; i1 < m_N1; i1++) {
598            for (dim_t i0 = 0; i0 < m_N0; i0++) {
599                double* point = arg.getSampleDataRW(i0+m_N0*i1);
600                point[0] = xdx.first+i0*xdx.second;
601                point[1] = ydy.first+i1*ydy.second;
602            }
603        }
604    }
605    
606    //protected
607    void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
608  {  {
     escript::Data& in = *const_cast<escript::Data*>(&cIn);  
609      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
610      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
611      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 302  void Rectangle::setToGradient(escript::D Line 627  void Rectangle::setToGradient(escript::D
627      const double cy7 = 1./h1;      const double cy7 = 1./h1;
628    
629      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
630          /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */          out.requireWrite();
631  #pragma omp parallel for  #pragma omp parallel for
632          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
633              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
634                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
635                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
636                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
637                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
638                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
639                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
640                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
# Line 323  void Rectangle::setToGradient(escript::D Line 648  void Rectangle::setToGradient(escript::D
648                  } /* end of component loop i */                  } /* end of component loop i */
649              } /* end of k0 loop */              } /* end of k0 loop */
650          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */  
651      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
652          /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */          out.requireWrite();
653  #pragma omp parallel for  #pragma omp parallel for
654          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
655              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
656                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
657                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
658                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
659                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
660                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
661                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
662                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
# Line 340  void Rectangle::setToGradient(escript::D Line 664  void Rectangle::setToGradient(escript::D
664                  } /* end of component loop i */                  } /* end of component loop i */
665              } /* end of k0 loop */              } /* end of k0 loop */
666          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */  
667      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
668            out.requireWrite();
669  #pragma omp parallel  #pragma omp parallel
670          {          {
             /*** GENERATOR SNIP_GRAD_FACES TOP */  
671              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
672  #pragma omp for nowait  #pragma omp for nowait
673                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
674                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
675                      const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
676                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
677                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
678                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
679                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
680                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
# Line 364  void Rectangle::setToGradient(escript::D Line 687  void Rectangle::setToGradient(escript::D
687              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
688  #pragma omp for nowait  #pragma omp for nowait
689                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
690                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
691                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
692                      const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
693                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
694                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
695                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
696                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
# Line 380  void Rectangle::setToGradient(escript::D Line 703  void Rectangle::setToGradient(escript::D
703              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
704  #pragma omp for nowait  #pragma omp for nowait
705                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
706                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
707                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
708                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
709                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
710                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
711                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
712                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
# Line 396  void Rectangle::setToGradient(escript::D Line 719  void Rectangle::setToGradient(escript::D
719              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
720  #pragma omp for nowait  #pragma omp for nowait
721                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
722                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
723                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
724                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
725                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
726                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
727                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
728                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
# Line 409  void Rectangle::setToGradient(escript::D Line 732  void Rectangle::setToGradient(escript::D
732                      } /* end of component loop i */                      } /* end of component loop i */
733                  } /* end of k0 loop */                  } /* end of k0 loop */
734              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_GRAD_FACES BOTTOM */  
735          } // end of parallel section          } // end of parallel section
736      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
737            out.requireWrite();
738  #pragma omp parallel  #pragma omp parallel
739          {          {
             /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */  
740              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
741  #pragma omp for nowait  #pragma omp for nowait
742                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
743                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
744                      const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
745                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
746                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
747                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
748                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
749                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
# Line 432  void Rectangle::setToGradient(escript::D Line 754  void Rectangle::setToGradient(escript::D
754              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
755  #pragma omp for nowait  #pragma omp for nowait
756                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
757                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
758                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
759                      const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
760                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
761                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
762                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
763                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
# Line 446  void Rectangle::setToGradient(escript::D Line 768  void Rectangle::setToGradient(escript::D
768              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
769  #pragma omp for nowait  #pragma omp for nowait
770                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
771                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
772                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
773                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
774                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
775                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
776                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
777                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
# Line 460  void Rectangle::setToGradient(escript::D Line 782  void Rectangle::setToGradient(escript::D
782              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
783  #pragma omp for nowait  #pragma omp for nowait
784                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
785                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
786                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
787                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
788                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
789                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
790                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
791                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
# Line 471  void Rectangle::setToGradient(escript::D Line 793  void Rectangle::setToGradient(escript::D
793                      } /* end of component loop i */                      } /* end of component loop i */
794                  } /* end of k0 loop */                  } /* end of k0 loop */
795              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */  
796          } // end of parallel section          } // end of parallel section
     } else {  
         stringstream msg;  
         msg << "setToGradient() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
797      }      }
798  }  }
799    
800  void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  //protected
801    void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
802  {  {
803      escript::Data& in = *const_cast<escript::Data*>(&arg);      const dim_t numComp = arg.getDataPointSize();
     const dim_t numComp = in.getDataPointSize();  
804      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
805      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
806        const index_t left = (m_offset0==0 ? 0 : 1);
807        const index_t bottom = (m_offset1==0 ? 0 : 1);
808      if (arg.getFunctionSpace().getTypeCode() == Elements) {      if (arg.getFunctionSpace().getTypeCode() == Elements) {
809          const double w_0 = h0*h1/4.;          const double w = h0*h1/4.;
810  #pragma omp parallel  #pragma omp parallel
811          {          {
812              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
813  #pragma omp for nowait  #pragma omp for nowait
814              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
815                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
816                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
817                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
818                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
819                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
820                          const register double f_2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
821                          const register double f_3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
822                          int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;                          int_local[i]+=(f0+f1+f2+f3)*w;
823                      }  /* end of component loop i */                      }  /* end of component loop i */
824                  } /* end of k0 loop */                  } /* end of k0 loop */
825              } /* end of k1 loop */              } /* end of k1 loop */
# Line 511  void Rectangle::setToIntegrals(vector<do Line 829  void Rectangle::setToIntegrals(vector<do
829                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
830          } // end of parallel section          } // end of parallel section
831      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
832          const double w_0 = h0*h1;          const double w = h0*h1;
833  #pragma omp parallel  #pragma omp parallel
834          {          {
835              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
836  #pragma omp for nowait  #pragma omp for nowait
837              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
838                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
839                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
840                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
841                          int_local[i]+=f[i]*w_0;                          int_local[i]+=f[i]*w;
842                      }  /* end of component loop i */                      }  /* end of component loop i */
843                  } /* end of k0 loop */                  } /* end of k0 loop */
844              } /* end of k1 loop */              } /* end of k1 loop */
# Line 530  void Rectangle::setToIntegrals(vector<do Line 848  void Rectangle::setToIntegrals(vector<do
848                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
849          } // end of parallel section          } // end of parallel section
850      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
851          const double w_0 = h0/2.;          const double w0 = h0/2.;
852          const double w_1 = h1/2.;          const double w1 = h1/2.;
853  #pragma omp parallel  #pragma omp parallel
854          {          {
855              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
856              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
857  #pragma omp for nowait  #pragma omp for nowait
858                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
859                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
860                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
861                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
862                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
863                          int_local[i]+=(f_0+f_1)*w_1;                          int_local[i]+=(f0+f1)*w1;
864                      }  /* end of component loop i */                      }  /* end of component loop i */
865                  } /* end of k1 loop */                  } /* end of k1 loop */
866              }              }
867    
868              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -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[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+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 double f0 = f[INDEX2(i,0,numComp)];
874                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const 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[2] > -1) {              if (m_faceOffset[2] > -1) {
881  #pragma omp for nowait  #pragma omp for nowait
882                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
883                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
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 double f0 = f[INDEX2(i,0,numComp)];
886                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
887                          int_local[i]+=(f_0+f_1)*w_0;                          int_local[i]+=(f0+f1)*w0;
888                      }  /* end of component loop i */                      }  /* end of component loop i */
889                  } /* end of k0 loop */                  } /* end of k0 loop */
890              }              }
891    
892              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -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[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+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 double f0 = f[INDEX2(i,0,numComp)];
898                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const 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              }              }
# Line 593  void Rectangle::setToIntegrals(vector<do Line 911  void Rectangle::setToIntegrals(vector<do
911              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
912              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
913  #pragma omp for nowait  #pragma omp for nowait
914                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
915                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
916                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
917                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
918                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 603  void Rectangle::setToIntegrals(vector<do Line 921  void Rectangle::setToIntegrals(vector<do
921    
922              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
923  #pragma omp for nowait  #pragma omp for nowait
924                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
925                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
926                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
927                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
928                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 613  void Rectangle::setToIntegrals(vector<do Line 931  void Rectangle::setToIntegrals(vector<do
931    
932              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
933  #pragma omp for nowait  #pragma omp for nowait
934                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
935                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
936                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
937                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
938                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 623  void Rectangle::setToIntegrals(vector<do Line 941  void Rectangle::setToIntegrals(vector<do
941    
942              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
943  #pragma omp for nowait  #pragma omp for nowait
944                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
945                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
946                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
947                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
948                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 635  void Rectangle::setToIntegrals(vector<do Line 953  void Rectangle::setToIntegrals(vector<do
953              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
954                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
955          } // 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;  
         }  
956      }      }
957  }  }
958    
# Line 1008  void Rectangle::createPattern() Line 1111  void Rectangle::createPattern()
1111      // The rest is assigned in the loop further down      // The rest is assigned in the loop further down
1112      m_dofMap.assign(getNumNodes(), 0);      m_dofMap.assign(getNumNodes(), 0);
1113  #pragma omp parallel for  #pragma omp parallel for
1114      for (index_t i=bottom; i<m_N1; i++) {      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1115          for (index_t j=left; j<m_N0; j++) {          for (index_t j=left; j<left+nDOF0; j++) {
1116              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
1117          }          }
1118      }      }
# Line 1170  void Rectangle::createPattern() Line 1273  void Rectangle::createPattern()
1273      Paso_Pattern_free(rowPattern);      Paso_Pattern_free(rowPattern);
1274  }  }
1275    
1276    //private
1277    void Rectangle::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
1278             const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
1279             bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
1280    {
1281        IndexVector rowIndex;
1282        rowIndex.push_back(m_dofMap[firstNode]);
1283        rowIndex.push_back(m_dofMap[firstNode+1]);
1284        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
1285        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
1286        if (addF) {
1287            double *F_p=F.getSampleDataRW(0);
1288            for (index_t i=0; i<rowIndex.size(); i++) {
1289                if (rowIndex[i]<getNumDOF()) {
1290                    for (index_t eq=0; eq<nEq; eq++) {
1291                        F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
1292                    }
1293                }
1294            }
1295        }
1296        if (addS) {
1297            addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
1298        }
1299    }
1300    
1301  //protected  //protected
1302  void Rectangle::interpolateNodesOnElements(escript::Data& out,  void Rectangle::interpolateNodesOnElements(escript::Data& out,
1303                                          escript::Data& in, bool reduced) const                                          escript::Data& in, bool reduced) const
1304  {  {
1305      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1306      if (reduced) {      if (reduced) {
1307          /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */          out.requireWrite();
1308          const double c0 = .25;          const double c0 = .25;
1309  #pragma omp parallel for  #pragma omp parallel for
1310          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1311              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
1312                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1313                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1314                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1315                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1316                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1317                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1318                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1319                  } /* end of component loop i */                  } /* end of component loop i */
1320              } /* end of k0 loop */              } /* end of k0 loop */
1321          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */  
1322      } else {      } else {
1323          /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */          out.requireWrite();
1324          const double c0 = .16666666666666666667;          const double c0 = .16666666666666666667;
1325          const double c1 = .044658198738520451079;          const double c1 = .044658198738520451079;
1326          const double c2 = .62200846792814621559;          const double c2 = .62200846792814621559;
1327  #pragma omp parallel for  #pragma omp parallel for
1328          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1329              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
1330                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1331                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1332                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1333                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1334                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1335                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1336                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);
# Line 1213  void Rectangle::interpolateNodesOnElemen Line 1340  void Rectangle::interpolateNodesOnElemen
1340                  } /* end of component loop i */                  } /* end of component loop i */
1341              } /* end of k0 loop */              } /* end of k0 loop */
1342          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */  
1343      }      }
1344  }  }
1345    
# Line 1223  void Rectangle::interpolateNodesOnFaces( Line 1349  void Rectangle::interpolateNodesOnFaces(
1349  {  {
1350      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1351      if (reduced) {      if (reduced) {
1352            out.requireWrite();
1353          const double c0 = .5;          const double c0 = .5;
1354  #pragma omp parallel  #pragma omp parallel
1355          {          {
             /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */  
1356              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1357  #pragma omp for nowait  #pragma omp for nowait
1358                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1359                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1360                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1361                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1362                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1363                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
# Line 1241  void Rectangle::interpolateNodesOnFaces( Line 1367  void Rectangle::interpolateNodesOnFaces(
1367              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1368  #pragma omp for nowait  #pragma omp for nowait
1369                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1370                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1371                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1372                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1373                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1374                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
# Line 1252  void Rectangle::interpolateNodesOnFaces( Line 1378  void Rectangle::interpolateNodesOnFaces(
1378              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1379  #pragma omp for nowait  #pragma omp for nowait
1380                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1381                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1382                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1383                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1384                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1385                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
# Line 1263  void Rectangle::interpolateNodesOnFaces( Line 1389  void Rectangle::interpolateNodesOnFaces(
1389              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1390  #pragma omp for nowait  #pragma omp for nowait
1391                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1392                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1393                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1394                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1395                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1396                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1397                      } /* end of component loop i */                      } /* end of component loop i */
1398                  } /* end of k0 loop */                  } /* end of k0 loop */
1399              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */  
1400          } // end of parallel section          } // end of parallel section
1401      } else {      } else {
1402            out.requireWrite();
1403          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1404          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1405  #pragma omp parallel  #pragma omp parallel
1406          {          {
             /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */  
1407              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1408  #pragma omp for nowait  #pragma omp for nowait
1409                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1410                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1411                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1412                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1413                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1414                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;
# Line 1294  void Rectangle::interpolateNodesOnFaces( Line 1419  void Rectangle::interpolateNodesOnFaces(
1419              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1420  #pragma omp for nowait  #pragma omp for nowait
1421                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1422                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1423                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1424                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1425                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1426                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;
# Line 1306  void Rectangle::interpolateNodesOnFaces( Line 1431  void Rectangle::interpolateNodesOnFaces(
1431              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1432  #pragma omp for nowait  #pragma omp for nowait
1433                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1434                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1435                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1436                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1437                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1438                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;
# Line 1318  void Rectangle::interpolateNodesOnFaces( Line 1443  void Rectangle::interpolateNodesOnFaces(
1443              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1444  #pragma omp for nowait  #pragma omp for nowait
1445                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1446                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1447                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1448                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1449                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1450                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;
# Line 1327  void Rectangle::interpolateNodesOnFaces( Line 1452  void Rectangle::interpolateNodesOnFaces(
1452                      } /* end of component loop i */                      } /* end of component loop i */
1453                  } /* end of k0 loop */                  } /* end of k0 loop */
1454              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */  
1455          } // end of parallel section          } // end of parallel section
1456      }      }
1457  }  }
# Line 1336  void Rectangle::interpolateNodesOnFaces( Line 1460  void Rectangle::interpolateNodesOnFaces(
1460  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1461          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1462          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
1463          const escript::Data& X, const escript::Data& Y,          const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
1464  {  {
1465      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1466      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
     /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */  
1467      const double w0 = -0.1555021169820365539*h1/h0;      const double w0 = -0.1555021169820365539*h1/h0;
1468      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
1469        const double w2 = -0.15550211698203655390;
1470        const double w3 = 0.041666666666666666667*h0/h1;
1471        const double w4 = 0.15550211698203655390;
1472        const double w5 = -0.041666666666666666667;
1473        const double w6 = -0.01116454968463011277*h1/h0;
1474        const double w7 = 0.011164549684630112770;
1475        const double w8 = -0.011164549684630112770;
1476        const double w9 = -0.041666666666666666667*h1/h0;
1477      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.041666666666666666667*h0/h1;
1478      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = 0.1555021169820365539*h1/h0;
1479      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = 0.1555021169820365539*h0/h1;
# Line 1353  void Rectangle::assemblePDESingle(Paso_S Line 1483  void Rectangle::assemblePDESingle(Paso_S
1483      const double w16 = -0.01116454968463011277*h0/h1;      const double w16 = -0.01116454968463011277*h0/h1;
1484      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*h0/h1;
1485      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*h1/h0;
1486      const double w19 = 0.25000000000000000000;      const double w19 = 0.25;
1487      const double w2 = -0.15550211698203655390;      const double w20 = -0.25;
     const double w20 = -0.25000000000000000000;  
1488      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*h0/h1;
1489      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*h1/h0;
1490      const double w23 = -0.16666666666666666667*h0/h1;      const double w23 = -0.16666666666666666667*h0/h1;
# Line 1365  void Rectangle::assemblePDESingle(Paso_S Line 1494  void Rectangle::assemblePDESingle(Paso_S
1494      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = -0.33333333333333333333*h0/h1;
1495      const double w28 = -0.032861463941450536761*h1;      const double w28 = -0.032861463941450536761*h1;
1496      const double w29 = -0.032861463941450536761*h0;      const double w29 = -0.032861463941450536761*h0;
     const double w3 = 0.041666666666666666667*h0/h1;  
1497      const double w30 = -0.12264065304058601714*h1;      const double w30 = -0.12264065304058601714*h1;
1498      const double w31 = -0.0023593469594139828636*h1;      const double w31 = -0.0023593469594139828636*h1;
1499      const double w32 = -0.008805202725216129906*h0;      const double w32 = -0.008805202725216129906*h0;
# Line 1376  void Rectangle::assemblePDESingle(Paso_S Line 1504  void Rectangle::assemblePDESingle(Paso_S
1504      const double w37 = 0.0023593469594139828636*h1;      const double w37 = 0.0023593469594139828636*h1;
1505      const double w38 = 0.12264065304058601714*h1;      const double w38 = 0.12264065304058601714*h1;
1506      const double w39 = 0.032861463941450536761*h0;      const double w39 = 0.032861463941450536761*h0;
     const double w4 = 0.15550211698203655390;  
1507      const double w40 = -0.12264065304058601714*h0;      const double w40 = -0.12264065304058601714*h0;
1508      const double w41 = -0.0023593469594139828636*h0;      const double w41 = -0.0023593469594139828636*h0;
1509      const double w42 = 0.0023593469594139828636*h0;      const double w42 = 0.0023593469594139828636*h0;
# Line 1387  void Rectangle::assemblePDESingle(Paso_S Line 1514  void Rectangle::assemblePDESingle(Paso_S
1514      const double w47 = 0.16666666666666666667*h1;      const double w47 = 0.16666666666666666667*h1;
1515      const double w48 = 0.083333333333333333333*h0;      const double w48 = 0.083333333333333333333*h0;
1516      const double w49 = -0.16666666666666666667*h0;      const double w49 = -0.16666666666666666667*h0;
     const double w5 = -0.041666666666666666667;  
1517      const double w50 = 0.16666666666666666667*h0;      const double w50 = 0.16666666666666666667*h0;
1518      const double w51 = -0.083333333333333333333*h1;      const double w51 = -0.083333333333333333333*h1;
1519      const double w52 = 0.025917019497006092316*h0*h1;      const double w52 = 0.025917019497006092316*h0*h1;
# Line 1398  void Rectangle::assemblePDESingle(Paso_S Line 1524  void Rectangle::assemblePDESingle(Paso_S
1524      const double w57 = 0.055555555555555555556*h0*h1;      const double w57 = 0.055555555555555555556*h0*h1;
1525      const double w58 = 0.027777777777777777778*h0*h1;      const double w58 = 0.027777777777777777778*h0*h1;
1526      const double w59 = 0.11111111111111111111*h0*h1;      const double w59 = 0.11111111111111111111*h0*h1;
     const double w6 = -0.01116454968463011277*h1/h0;  
1527      const double w60 = -0.19716878364870322056*h1;      const double w60 = -0.19716878364870322056*h1;
1528      const double w61 = -0.19716878364870322056*h0;      const double w61 = -0.19716878364870322056*h0;
1529      const double w62 = -0.052831216351296779436*h0;      const double w62 = -0.052831216351296779436*h0;
# Line 1409  void Rectangle::assemblePDESingle(Paso_S Line 1534  void Rectangle::assemblePDESingle(Paso_S
1534      const double w67 = 0.052831216351296779436*h0;      const double w67 = 0.052831216351296779436*h0;
1535      const double w68 = -0.5*h1;      const double w68 = -0.5*h1;
1536      const double w69 = -0.5*h0;      const double w69 = -0.5*h0;
     const double w7 = 0.011164549684630112770;  
1537      const double w70 = 0.5*h1;      const double w70 = 0.5*h1;
1538      const double w71 = 0.5*h0;      const double w71 = 0.5*h0;
1539      const double w72 = 0.1555021169820365539*h0*h1;      const double w72 = 0.1555021169820365539*h0*h1;
1540      const double w73 = 0.041666666666666666667*h0*h1;      const double w73 = 0.041666666666666666667*h0*h1;
1541      const double w74 = 0.01116454968463011277*h0*h1;      const double w74 = 0.01116454968463011277*h0*h1;
1542      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 */  
1543    
1544      rhs.requireWrite();      rhs.requireWrite();
1545  #pragma omp parallel  #pragma omp parallel
# Line 1432  void Rectangle::assemblePDESingle(Paso_S Line 1553  void Rectangle::assemblePDESingle(Paso_S
1553                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1554                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1555                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
                     /*** GENERATOR SNIP_PDE_SINGLE TOP */  
1556                      ///////////////                      ///////////////
1557                      // process A //                      // process A //
1558                      ///////////////                      ///////////////
# Line 1440  void Rectangle::assemblePDESingle(Paso_S Line 1560  void Rectangle::assemblePDESingle(Paso_S
1560                          add_EM_S=true;                          add_EM_S=true;
1561                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1562                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
1563                              const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];                              const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1564                              const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];                              const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1565                              const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];                              const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1566                              const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];                              const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1567                              const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];                              const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1568                              const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];                              const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1569                              const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];                              const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1570                              const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];                              const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1571                              const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];                              const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1572                              const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];                              const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1573                              const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];                              const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1574                              const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];                              const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1575                              const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];                              const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1576                              const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];                              const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1577                              const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];                              const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1578                              const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];                              const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1579                              const register double tmp4_0 = A_10_1 + A_10_2;                              const double tmp0_0 = A_01_0 + A_01_3;
1580                              const register double tmp12_0 = A_11_0 + A_11_2;                              const double tmp1_0 = A_00_0 + A_00_1;
1581                              const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                              const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1582                              const register double tmp10_0 = A_01_3 + A_10_3;                              const double tmp3_0 = A_00_2 + A_00_3;
1583                              const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                              const double tmp4_0 = A_10_1 + A_10_2;
1584                              const register double tmp0_0 = A_01_0 + A_01_3;                              const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1585                              const register double tmp13_0 = A_01_2 + A_10_1;                              const double tmp6_0 = A_01_3 + A_10_0;
1586                              const register double tmp3_0 = A_00_2 + A_00_3;                              const double tmp7_0 = A_01_0 + A_10_3;
1587                              const register double tmp11_0 = A_11_1 + A_11_3;                              const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1588                              const register double tmp18_0 = A_01_1 + A_10_1;                              const double tmp9_0 = A_01_0 + A_10_0;
1589                              const register double tmp1_0 = A_00_0 + A_00_1;                              const double tmp12_0 = A_11_0 + A_11_2;
1590                              const register double tmp15_0 = A_01_1 + A_10_2;                              const double tmp10_0 = A_01_3 + A_10_3;
1591                              const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                              const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
1592                              const register double tmp16_0 = A_10_0 + A_10_3;                              const double tmp13_0 = A_01_2 + A_10_1;
1593                              const register double tmp6_0 = A_01_3 + A_10_0;                              const double tmp11_0 = A_11_1 + A_11_3;
1594                              const register double tmp17_0 = A_01_1 + A_01_2;                              const double tmp18_0 = A_01_1 + A_10_1;
1595                              const register double tmp9_0 = A_01_0 + A_10_0;                              const double tmp15_0 = A_01_1 + A_10_2;
1596                              const register double tmp7_0 = A_01_0 + A_10_3;                              const double tmp16_0 = A_10_0 + A_10_3;
1597                              const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                              const double tmp17_0 = A_01_1 + A_01_2;
1598                              const register double tmp19_0 = A_01_2 + A_10_2;                              const double tmp19_0 = A_01_2 + A_10_2;
1599                              const register double tmp14_1 = A_10_0*w8;                              const double tmp0_1 = A_10_3*w8;
1600                              const register double tmp23_1 = tmp3_0*w14;                              const double tmp1_1 = tmp0_0*w1;
1601                              const register double tmp35_1 = A_01_0*w8;                              const double tmp2_1 = A_01_1*w4;
1602                              const register double tmp54_1 = tmp13_0*w8;                              const double tmp3_1 = tmp1_0*w0;
1603                              const register double tmp20_1 = tmp9_0*w4;                              const double tmp4_1 = A_01_2*w7;
1604                              const register double tmp25_1 = tmp12_0*w12;                              const double tmp5_1 = tmp2_0*w3;
1605                              const register double tmp2_1 = A_01_1*w4;                              const double tmp6_1 = tmp3_0*w6;
1606                              const register double tmp44_1 = tmp7_0*w7;                              const double tmp7_1 = A_10_0*w2;
1607                              const register double tmp26_1 = tmp10_0*w4;                              const double tmp8_1 = tmp4_0*w5;
1608                              const register double tmp52_1 = tmp18_0*w8;                              const double tmp9_1 = tmp2_0*w10;
1609                              const register double tmp48_1 = A_10_1*w7;                              const double tmp14_1 = A_10_0*w8;
1610                              const register double tmp46_1 = A_01_3*w8;                              const double tmp23_1 = tmp3_0*w14;
1611                              const register double tmp50_1 = A_01_0*w2;                              const double tmp35_1 = A_01_0*w8;
1612                              const register double tmp8_1 = tmp4_0*w5;                              const double tmp54_1 = tmp13_0*w8;
1613                              const register double tmp56_1 = tmp19_0*w8;                              const double tmp20_1 = tmp9_0*w4;
1614                              const register double tmp9_1 = tmp2_0*w10;                              const double tmp25_1 = tmp12_0*w12;
1615                              const register double tmp19_1 = A_10_3*w2;                              const double tmp44_1 = tmp7_0*w7;
1616                              const register double tmp47_1 = A_10_2*w4;                              const double tmp26_1 = tmp10_0*w4;
1617                              const register double tmp16_1 = tmp3_0*w0;                              const double tmp52_1 = tmp18_0*w8;
1618                              const register double tmp18_1 = tmp1_0*w6;                              const double tmp48_1 = A_10_1*w7;
1619                              const register double tmp31_1 = tmp11_0*w12;                              const double tmp46_1 = A_01_3*w8;
1620                              const register double tmp55_1 = tmp15_0*w2;                              const double tmp50_1 = A_01_0*w2;
1621                              const register double tmp39_1 = A_10_2*w7;                              const double tmp56_1 = tmp19_0*w8;
1622                              const register double tmp11_1 = tmp6_0*w7;                              const double tmp19_1 = A_10_3*w2;
1623                              const register double tmp40_1 = tmp11_0*w17;                              const double tmp47_1 = A_10_2*w4;
1624                              const register double tmp34_1 = tmp15_0*w8;                              const double tmp16_1 = tmp3_0*w0;
1625                              const register double tmp33_1 = tmp14_0*w5;                              const double tmp18_1 = tmp1_0*w6;
1626                              const register double tmp24_1 = tmp11_0*w13;                              const double tmp31_1 = tmp11_0*w12;
1627                              const register double tmp3_1 = tmp1_0*w0;                              const double tmp55_1 = tmp15_0*w2;
1628                              const register double tmp5_1 = tmp2_0*w3;                              const double tmp39_1 = A_10_2*w7;
1629                              const register double tmp43_1 = tmp17_0*w5;                              const double tmp11_1 = tmp6_0*w7;
1630                              const register double tmp15_1 = A_01_2*w4;                              const double tmp40_1 = tmp11_0*w17;
1631                              const register double tmp53_1 = tmp19_0*w2;                              const double tmp34_1 = tmp15_0*w8;
1632                              const register double tmp27_1 = tmp3_0*w11;                              const double tmp33_1 = tmp14_0*w5;
1633                              const register double tmp32_1 = tmp13_0*w2;                              const double tmp24_1 = tmp11_0*w13;
1634                              const register double tmp10_1 = tmp5_0*w9;                              const double tmp43_1 = tmp17_0*w5;
1635                              const register double tmp37_1 = A_10_1*w4;                              const double tmp15_1 = A_01_2*w4;
1636                              const register double tmp38_1 = tmp5_0*w15;                              const double tmp53_1 = tmp19_0*w2;
1637                              const register double tmp17_1 = A_01_1*w7;                              const double tmp27_1 = tmp3_0*w11;
1638                              const register double tmp12_1 = tmp7_0*w4;                              const double tmp32_1 = tmp13_0*w2;
1639                              const register double tmp22_1 = tmp10_0*w7;                              const double tmp10_1 = tmp5_0*w9;
1640                              const register double tmp57_1 = tmp18_0*w2;                              const double tmp37_1 = A_10_1*w4;
1641                              const register double tmp28_1 = tmp9_0*w7;                              const double tmp38_1 = tmp5_0*w15;
1642                              const register double tmp29_1 = tmp1_0*w14;                              const double tmp17_1 = A_01_1*w7;
1643                              const register double tmp51_1 = tmp11_0*w16;                              const double tmp12_1 = tmp7_0*w4;
1644                              const register double tmp42_1 = tmp12_0*w16;                              const double tmp22_1 = tmp10_0*w7;
1645                              const register double tmp49_1 = tmp12_0*w17;                              const double tmp57_1 = tmp18_0*w2;
1646                              const register double tmp21_1 = tmp1_0*w11;                              const double tmp28_1 = tmp9_0*w7;
1647                              const register double tmp1_1 = tmp0_0*w1;                              const double tmp29_1 = tmp1_0*w14;
1648                              const register double tmp45_1 = tmp6_0*w4;                              const double tmp51_1 = tmp11_0*w16;
1649                              const register double tmp7_1 = A_10_0*w2;                              const double tmp42_1 = tmp12_0*w16;
1650                              const register double tmp6_1 = tmp3_0*w6;                              const double tmp49_1 = tmp12_0*w17;
1651                              const register double tmp13_1 = tmp8_0*w1;                              const double tmp21_1 = tmp1_0*w11;
1652                              const register double tmp36_1 = tmp16_0*w1;                              const double tmp45_1 = tmp6_0*w4;
1653                              const register double tmp41_1 = A_01_3*w2;                              const double tmp13_1 = tmp8_0*w1;
1654                              const register double tmp30_1 = tmp12_0*w13;                              const double tmp36_1 = tmp16_0*w1;
1655                              const register double tmp4_1 = A_01_2*w7;                              const double tmp41_1 = A_01_3*w2;
1656                              const register double tmp0_1 = A_10_3*w8;                              const double tmp30_1 = tmp12_0*w13;
                             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;  
1657                              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;
1658                              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;
1659                                EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
1660                              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;
1661                              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;
1662                                EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
1663                              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;
1664                                EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
1665                              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;
1666                              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;  
1667                              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;
1668                              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;
1669                              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;
1670                              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;
1671                          } 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;
1672                              const register double A_00 = A_p[INDEX2(0,0,2)];                              EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1673                              const register double A_01 = A_p[INDEX2(0,1,2)];                          } else { // constant data
1674                              const register double A_10 = A_p[INDEX2(1,0,2)];                              const double A_00 = A_p[INDEX2(0,0,2)];
1675                              const register double A_11 = A_p[INDEX2(1,1,2)];                              const double A_10 = A_p[INDEX2(1,0,2)];
1676                              const register double tmp0_0 = A_01 + A_10;                              const double A_01 = A_p[INDEX2(0,1,2)];
1677                              const register double tmp0_1 = A_00*w18;                              const double A_11 = A_p[INDEX2(1,1,2)];
1678                              const register double tmp10_1 = A_01*w20;                              const double tmp0_0 = A_01 + A_10;
1679                              const register double tmp12_1 = A_00*w26;                              const double tmp0_1 = A_00*w18;
1680                              const register double tmp4_1 = A_00*w22;                              const double tmp1_1 = A_01*w19;
1681                              const register double tmp8_1 = A_00*w24;                              const double tmp2_1 = A_10*w20;
1682                              const register double tmp13_1 = A_10*w19;                              const double tmp3_1 = A_11*w21;
1683                              const register double tmp9_1 = tmp0_0*w20;                              const double tmp4_1 = A_00*w22;
1684                              const register double tmp3_1 = A_11*w21;                              const double tmp5_1 = tmp0_0*w19;
1685                              const register double tmp11_1 = A_11*w27;                              const double tmp6_1 = A_11*w23;
1686                              const register double tmp1_1 = A_01*w19;                              const double tmp7_1 = A_11*w25;
1687                              const register double tmp6_1 = A_11*w23;                              const double tmp8_1 = A_00*w24;
1688                              const register double tmp7_1 = A_11*w25;                              const double tmp9_1 = tmp0_0*w20;
1689                              const register double tmp2_1 = A_10*w20;                              const double tmp10_1 = A_01*w20;
1690                              const register double tmp5_1 = tmp0_0*w19;                              const double tmp11_1 = A_11*w27;
1691                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              const double tmp12_1 = A_00*w26;
1692                              EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;                              const double tmp13_1 = A_10*w19;
                             EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
1693                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1694                              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;
1695                                EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1696                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1697                              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;
1698                                EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1699                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1700                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1701                              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;
1702                              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;  
1703                              EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1704                              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;
1705                              EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1706                              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;
1707                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1708                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1709                          }                          }
1710                      }                      }
1711                      ///////////////                      ///////////////
# Line 1595  void Rectangle::assemblePDESingle(Paso_S Line 1715  void Rectangle::assemblePDESingle(Paso_S
1715                          add_EM_S=true;                          add_EM_S=true;
1716                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1717                          if (B.actsExpanded()) {                          if (B.actsExpanded()) {
1718                              const register double B_0_0 = B_p[INDEX2(0,0,2)];                              const double B_0_0 = B_p[INDEX2(0,0,2)];
1719                              const register double B_1_0 = B_p[INDEX2(1,0,2)];                              const double B_1_0 = B_p[INDEX2(1,0,2)];
1720                              const register double B_0_1 = B_p[INDEX2(0,1,2)];                              const double B_0_1 = B_p[INDEX2(0,1,2)];
1721                              const register double B_1_1 = B_p[INDEX2(1,1,2)];                              const double B_1_1 = B_p[INDEX2(1,1,2)];
1722                              const register double B_0_2 = B_p[INDEX2(0,2,2)];                              const double B_0_2 = B_p[INDEX2(0,2,2)];
1723                              const register double B_1_2 = B_p[INDEX2(1,2,2)];                              const double B_1_2 = B_p[INDEX2(1,2,2)];
1724                              const register double B_0_3 = B_p[INDEX2(0,3,2)];                              const double B_0_3 = B_p[INDEX2(0,3,2)];
1725                              const register double B_1_3 = B_p[INDEX2(1,3,2)];                              const double B_1_3 = B_p[INDEX2(1,3,2)];
1726                              const register double tmp3_0 = B_0_0 + B_0_2;                              const double tmp0_0 = B_1_0 + B_1_1;
1727                              const register double tmp1_0 = B_1_2 + B_1_3;                              const double tmp1_0 = B_1_2 + B_1_3;
1728                              const register double tmp2_0 = B_0_1 + B_0_3;                              const double tmp2_0 = B_0_1 + B_0_3;
1729                              const register double tmp0_0 = B_1_0 + B_1_1;                              const double tmp3_0 = B_0_0 + B_0_2;
1730                              const register double tmp63_1 = B_1_1*w42;                              const double tmp63_1 = B_1_1*w42;
1731                              const register double tmp79_1 = B_1_1*w40;                              const double tmp79_1 = B_1_1*w40;
1732                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
1733                              const register double tmp8_1 = tmp0_0*w32;                              const double tmp8_1 = tmp0_0*w32;
1734                              const register double tmp71_1 = B_0_1*w34;                              const double tmp71_1 = B_0_1*w34;
1735                              const register double tmp19_1 = B_0_3*w31;                              const double tmp19_1 = B_0_3*w31;
1736                              const register double tmp15_1 = B_0_3*w34;                              const double tmp15_1 = B_0_3*w34;
1737                              const register double tmp9_1 = tmp3_0*w34;                              const double tmp9_1 = tmp3_0*w34;
1738                              const register double tmp35_1 = B_1_0*w36;                              const double tmp35_1 = B_1_0*w36;
1739                              const register double tmp66_1 = B_0_3*w28;                              const double tmp66_1 = B_0_3*w28;
1740                              const register double tmp28_1 = B_1_0*w42;                              const double tmp28_1 = B_1_0*w42;
1741                              const register double tmp22_1 = B_1_0*w40;                              const double tmp22_1 = B_1_0*w40;
1742                              const register double tmp16_1 = B_1_2*w29;                              const double tmp16_1 = B_1_2*w29;
1743                              const register double tmp6_1 = tmp2_0*w35;                              const double tmp6_1 = tmp2_0*w35;
1744                              const register double tmp55_1 = B_1_3*w40;                              const double tmp55_1 = B_1_3*w40;
1745                              const register double tmp50_1 = B_1_3*w42;                              const double tmp50_1 = B_1_3*w42;
1746                              const register double tmp7_1 = tmp1_0*w29;                              const double tmp7_1 = tmp1_0*w29;
1747                              const register double tmp1_1 = tmp1_0*w32;                              const double tmp1_1 = tmp1_0*w32;
1748                              const register double tmp57_1 = B_0_3*w30;                              const double tmp57_1 = B_0_3*w30;
1749                              const register double tmp18_1 = B_1_1*w32;                              const double tmp18_1 = B_1_1*w32;
1750                              const register double tmp53_1 = B_1_0*w41;                              const double tmp53_1 = B_1_0*w41;
1751                              const register double tmp61_1 = B_1_3*w36;                              const double tmp61_1 = B_1_3*w36;
1752                              const register double tmp27_1 = B_0_3*w38;                              const double tmp27_1 = B_0_3*w38;
1753                              const register double tmp64_1 = B_0_2*w30;                              const double tmp64_1 = B_0_2*w30;
1754                              const register double tmp76_1 = B_0_1*w38;                              const double tmp76_1 = B_0_1*w38;
1755                              const register double tmp39_1 = tmp2_0*w34;                              const double tmp39_1 = tmp2_0*w34;
1756                              const register double tmp62_1 = B_0_1*w31;                              const double tmp62_1 = B_0_1*w31;
1757                              const register double tmp56_1 = B_0_0*w31;                              const double tmp56_1 = B_0_0*w31;
1758                              const register double tmp49_1 = B_1_1*w36;                              const double tmp49_1 = B_1_1*w36;
1759                              const register double tmp2_1 = B_0_2*w31;                              const double tmp2_1 = B_0_2*w31;
1760                              const register double tmp23_1 = B_0_2*w33;                              const double tmp23_1 = B_0_2*w33;
1761                              const register double tmp38_1 = B_1_1*w43;                              const double tmp38_1 = B_1_1*w43;
1762                              const register double tmp74_1 = B_1_2*w41;                              const double tmp74_1 = B_1_2*w41;
1763                              const register double tmp43_1 = B_1_1*w41;                              const double tmp43_1 = B_1_1*w41;
1764                              const register double tmp58_1 = B_0_2*w28;                              const double tmp58_1 = B_0_2*w28;
1765                              const register double tmp67_1 = B_0_0*w33;                              const double tmp67_1 = B_0_0*w33;
1766                              const register double tmp33_1 = tmp0_0*w39;                              const double tmp33_1 = tmp0_0*w39;
1767                              const register double tmp4_1 = B_0_0*w28;                              const double tmp4_1 = B_0_0*w28;
1768                              const register double tmp20_1 = B_0_0*w30;                              const double tmp20_1 = B_0_0*w30;
1769                              const register double tmp13_1 = B_0_2*w38;                              const double tmp13_1 = B_0_2*w38;
1770                              const register double tmp65_1 = B_1_2*w43;                              const double tmp65_1 = B_1_2*w43;
1771                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
1772                              const register double tmp41_1 = tmp3_0*w33;                              const double tmp41_1 = tmp3_0*w33;
1773                              const register double tmp73_1 = B_0_2*w37;                              const double tmp73_1 = B_0_2*w37;
1774                              const register double tmp69_1 = B_0_0*w38;                              const double tmp69_1 = B_0_0*w38;
1775                              const register double tmp48_1 = B_1_2*w39;                              const double tmp48_1 = B_1_2*w39;
1776                              const register double tmp59_1 = B_0_1*w33;                              const double tmp59_1 = B_0_1*w33;
1777                              const register double tmp17_1 = B_1_3*w41;                              const double tmp17_1 = B_1_3*w41;
1778                              const register double tmp5_1 = B_0_3*w33;                              const double tmp5_1 = B_0_3*w33;
1779                              const register double tmp3_1 = B_0_1*w30;                              const double tmp3_1 = B_0_1*w30;
1780                              const register double tmp21_1 = B_0_1*w28;                              const double tmp21_1 = B_0_1*w28;
1781                              const register double tmp42_1 = B_1_0*w29;                              const double tmp42_1 = B_1_0*w29;
1782                              const register double tmp54_1 = B_1_2*w32;                              const double tmp54_1 = B_1_2*w32;
1783                              const register double tmp60_1 = B_1_0*w39;                              const double tmp60_1 = B_1_0*w39;
1784                              const register double tmp32_1 = tmp1_0*w36;                              const double tmp32_1 = tmp1_0*w36;
1785                              const register double tmp10_1 = B_0_1*w37;                              const double tmp10_1 = B_0_1*w37;
1786                              const register double tmp14_1 = B_0_0*w35;                              const double tmp14_1 = B_0_0*w35;
1787                              const register double tmp29_1 = B_0_1*w35;                              const double tmp29_1 = B_0_1*w35;
1788                              const register double tmp26_1 = B_1_2*w36;                              const double tmp26_1 = B_1_2*w36;
1789                              const register double tmp30_1 = B_1_3*w43;                              const double tmp30_1 = B_1_3*w43;
1790                              const register double tmp70_1 = B_0_2*w35;                              const double tmp70_1 = B_0_2*w35;
1791                              const register double tmp34_1 = B_1_3*w39;                              const double tmp34_1 = B_1_3*w39;
1792                              const register double tmp51_1 = B_1_0*w43;                              const double tmp51_1 = B_1_0*w43;
1793                              const register double tmp31_1 = B_0_2*w34;                              const double tmp31_1 = B_0_2*w34;
1794                              const register double tmp45_1 = tmp3_0*w28;                              const double tmp45_1 = tmp3_0*w28;
1795                              const register double tmp11_1 = tmp1_0*w39;                              const double tmp11_1 = tmp1_0*w39;
1796                              const register double tmp52_1 = B_1_1*w29;                              const double tmp52_1 = B_1_1*w29;
1797                              const register double tmp44_1 = B_1_3*w32;                              const double tmp44_1 = B_1_3*w32;
1798                              const register double tmp25_1 = B_1_1*w39;                              const double tmp25_1 = B_1_1*w39;
1799                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
1800                              const register double tmp72_1 = B_1_3*w29;                              const double tmp72_1 = B_1_3*w29;
1801                              const register double tmp40_1 = tmp2_0*w28;                              const double tmp40_1 = tmp2_0*w28;
1802                              const register double tmp46_1 = B_1_2*w40;                              const double tmp46_1 = B_1_2*w40;
1803                              const register double tmp36_1 = B_1_2*w42;                              const double tmp36_1 = B_1_2*w42;
1804                              const register double tmp24_1 = B_0_0*w37;                              const double tmp24_1 = B_0_0*w37;
1805                              const register double tmp77_1 = B_0_3*w35;                              const double tmp77_1 = B_0_3*w35;
1806                              const register double tmp68_1 = B_0_3*w37;                              const double tmp68_1 = B_0_3*w37;
1807                              const register double tmp78_1 = B_0_0*w34;                              const double tmp78_1 = B_0_0*w34;
1808                              const register double tmp12_1 = tmp0_0*w36;                              const double tmp12_1 = tmp0_0*w36;
1809                              const register double tmp75_1 = B_1_0*w32;                              const 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;  
1810                              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;
1811                              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;
1812                                EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1813                              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;
1814                              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;
1815                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1816                              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;
1817                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1818                              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;
1819                              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;  
1820                              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;
1821                              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;
1822                              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;
1823                              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;
1824                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1825                              const register double B_0 = B_p[0];                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1826                              const register double B_1 = B_p[1];                          } else { // constant data
1827                              const register double tmp6_1 = B_1*w50;                              const double B_0 = B_p[0];
1828                              const register double tmp1_1 = B_1*w45;                              const double B_1 = B_p[1];
1829                              const register double tmp5_1 = B_1*w49;                              const double tmp0_1 = B_0*w44;
1830                              const register double tmp4_1 = B_1*w48;                              const double tmp1_1 = B_1*w45;
1831                              const register double tmp0_1 = B_0*w44;                              const double tmp2_1 = B_0*w46;
1832                              const register double tmp2_1 = B_0*w46;                              const double tmp3_1 = B_0*w47;
1833                              const register double tmp7_1 = B_0*w51;                              const double tmp4_1 = B_1*w48;
1834                              const register double tmp3_1 = B_0*w47;                              const double tmp5_1 = B_1*w49;
1835                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const double tmp6_1 = B_1*w50;
1836                              EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;                              const double tmp7_1 = B_0*w51;
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
1837                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1838                              EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1839                                EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
1840                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
1841                              EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1842                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
1843                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
1844                                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
1845                              EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;                              EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
1846                              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;  
1847                              EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;                              EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
1848                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1849                              EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;                              EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
1850                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
1851                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
1852                                EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
1853                          }                          }
1854                      }                      }
1855                      ///////////////                      ///////////////
# Line 1739  void Rectangle::assemblePDESingle(Paso_S Line 1859  void Rectangle::assemblePDESingle(Paso_S
1859                          add_EM_S=true;                          add_EM_S=true;
1860                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
1861                          if (C.actsExpanded()) {                          if (C.actsExpanded()) {
1862                              const register double C_0_0 = C_p[INDEX2(0,0,2)];                              const double C_0_0 = C_p[INDEX2(0,0,2)];
1863                              const register double C_1_0 = C_p[INDEX2(1,0,2)];                              const double C_1_0 = C_p[INDEX2(1,0,2)];
1864                              const register double C_0_1 = C_p[INDEX2(0,1,2)];                              const double C_0_1 = C_p[INDEX2(0,1,2)];
1865                              const register double C_1_1 = C_p[INDEX2(1,1,2)];                              const double C_1_1 = C_p[INDEX2(1,1,2)];
1866                              const register double C_0_2 = C_p[INDEX2(0,2,2)];                              const double C_0_2 = C_p[INDEX2(0,2,2)];
1867                              const register double C_1_2 = C_p[INDEX2(1,2,2)];                              const double C_1_2 = C_p[INDEX2(1,2,2)];
1868                              const register double C_0_3 = C_p[INDEX2(0,3,2)];                              const double C_0_3 = C_p[INDEX2(0,3,2)];
1869                              const register double C_1_3 = C_p[INDEX2(1,3,2)];                              const double C_1_3 = C_p[INDEX2(1,3,2)];
1870                              const register double tmp2_0 = C_0_1 + C_0_3;                              const double tmp0_0 = C_1_0 + C_1_1;
1871                              const register double tmp1_0 = C_1_2 + C_1_3;                              const double tmp1_0 = C_1_2 + C_1_3;
1872                              const register double tmp3_0 = C_0_0 + C_0_2;                              const double tmp2_0 = C_0_1 + C_0_3;
1873                              const register double tmp0_0 = C_1_0 + C_1_1;                              const double tmp3_0 = C_0_0 + C_0_2;
1874                              const register double tmp64_1 = C_0_2*w30;                              const double tmp64_1 = C_0_2*w30;
1875                              const register double tmp14_1 = C_0_2*w28;                              const double tmp14_1 = C_0_2*w28;
1876                              const register double tmp19_1 = C_0_3*w31;                              const double tmp19_1 = C_0_3*w31;
1877                              const register double tmp22_1 = C_1_0*w40;                              const double tmp22_1 = C_1_0*w40;
1878                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
1879                              const register double tmp29_1 = C_0_1*w35;                              const double tmp29_1 = C_0_1*w35;
1880                              const register double tmp73_1 = C_0_2*w37;                              const double tmp73_1 = C_0_2*w37;
1881                              const register double tmp74_1 = C_1_2*w41;                              const double tmp74_1 = C_1_2*w41;
1882                              const register double tmp52_1 = C_1_3*w39;                              const double tmp52_1 = C_1_3*w39;
1883                              const register double tmp25_1 = C_1_1*w39;                              const double tmp25_1 = C_1_1*w39;
1884                              const register double tmp62_1 = C_0_1*w31;                              const double tmp62_1 = C_0_1*w31;
1885                              const register double tmp79_1 = C_1_1*w40;                              const double tmp79_1 = C_1_1*w40;
1886                              const register double tmp43_1 = C_1_1*w36;                              const double tmp43_1 = C_1_1*w36;
1887                              const register double tmp27_1 = C_0_3*w38;                              const double tmp27_1 = C_0_3*w38;
1888                              const register double tmp28_1 = C_1_0*w42;                              const double tmp28_1 = C_1_0*w42;
1889                              const register double tmp63_1 = C_1_1*w42;                              const double tmp63_1 = C_1_1*w42;
1890                              const register double tmp59_1 = C_0_3*w34;                              const double tmp59_1 = C_0_3*w34;
1891                              const register double tmp72_1 = C_1_3*w29;                              const double tmp72_1 = C_1_3*w29;
1892                              const register double tmp40_1 = tmp2_0*w35;                              const double tmp40_1 = tmp2_0*w35;
1893                              const register double tmp13_1 = C_0_3*w30;                              const double tmp13_1 = C_0_3*w30;
1894                              const register double tmp51_1 = C_1_2*w40;                              const double tmp51_1 = C_1_2*w40;
1895                              const register double tmp54_1 = C_1_2*w42;                              const double tmp54_1 = C_1_2*w42;
1896                              const register double tmp12_1 = C_0_0*w31;                              const double tmp12_1 = C_0_0*w31;
1897                              const register double tmp2_1 = tmp1_0*w32;                              const double tmp2_1 = tmp1_0*w32;
1898                              const register double tmp68_1 = C_0_2*w31;                              const double tmp68_1 = C_0_2*w31;
1899                              const register double tmp75_1 = C_1_0*w32;                              const double tmp75_1 = C_1_0*w32;
1900                              const register double tmp49_1 = C_1_1*w41;                              const double tmp49_1 = C_1_1*w41;
1901                              const register double tmp4_1 = C_0_2*w35;                              const double tmp4_1 = C_0_2*w35;
1902                              const register double tmp66_1 = C_0_3*w28;                              const double tmp66_1 = C_0_3*w28;
1903                              const register double tmp56_1 = C_0_1*w37;                              const double tmp56_1 = C_0_1*w37;
1904                              const register double tmp5_1 = C_0_1*w34;                              const double tmp5_1 = C_0_1*w34;
1905                              const register double tmp38_1 = tmp2_0*w34;                              const double tmp38_1 = tmp2_0*w34;
1906                              const register double tmp76_1 = C_0_1*w38;                              const double tmp76_1 = C_0_1*w38;
1907                              const register double tmp21_1 = C_0_1*w28;                              const double tmp21_1 = C_0_1*w28;
1908                              const register double tmp69_1 = C_0_1*w30;                              const double tmp69_1 = C_0_1*w30;
1909                              const register double tmp53_1 = C_1_0*w36;                              const double tmp53_1 = C_1_0*w36;
1910                              const register double tmp42_1 = C_1_2*w39;                              const double tmp42_1 = C_1_2*w39;
1911                              const register double tmp32_1 = tmp1_0*w29;                              const double tmp32_1 = tmp1_0*w29;
1912                              const register double tmp45_1 = C_1_0*w43;                              const double tmp45_1 = C_1_0*w43;
1913                              const register double tmp33_1 = tmp0_0*w32;                              const double tmp33_1 = tmp0_0*w32;
1914                              const register double tmp35_1 = C_1_0*w41;                              const double tmp35_1 = C_1_0*w41;
1915                              const register double tmp26_1 = C_1_2*w36;                              const double tmp26_1 = C_1_2*w36;
1916                              const register double tmp67_1 = C_0_0*w33;                              const double tmp67_1 = C_0_0*w33;
1917                              const register double tmp31_1 = C_0_2*w34;                              const double tmp31_1 = C_0_2*w34;
1918                              const register double tmp20_1 = C_0_0*w30;                              const double tmp20_1 = C_0_0*w30;
1919                              const register double tmp70_1 = C_0_0*w28;                              const double tmp70_1 = C_0_0*w28;
1920                              const register double tmp8_1 = tmp0_0*w39;                              const double tmp8_1 = tmp0_0*w39;
1921                              const register double tmp30_1 = C_1_3*w43;                              const double tmp30_1 = C_1_3*w43;
1922                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
1923                              const register double tmp17_1 = C_1_3*w41;                              const double tmp17_1 = C_1_3*w41;
1924                              const register double tmp58_1 = C_0_0*w35;                              const double tmp58_1 = C_0_0*w35;
1925                              const register double tmp9_1 = tmp3_0*w33;                              const double tmp9_1 = tmp3_0*w33;
1926                              const register double tmp61_1 = C_1_3*w36;                              const double tmp61_1 = C_1_3*w36;
1927                              const register double tmp41_1 = tmp3_0*w34;                              const double tmp41_1 = tmp3_0*w34;
1928                              const register double tmp50_1 = C_1_3*w32;                              const double tmp50_1 = C_1_3*w32;
1929                              const register double tmp18_1 = C_1_1*w32;                              const double tmp18_1 = C_1_1*w32;
1930                              const register double tmp6_1 = tmp1_0*w36;                              const double tmp6_1 = tmp1_0*w36;
1931                              const register double tmp3_1 = C_0_0*w38;                              const double tmp3_1 = C_0_0*w38;
1932                              const register double tmp34_1 = C_1_1*w29;                              const double tmp34_1 = C_1_1*w29;
1933                              const register double tmp77_1 = C_0_3*w35;                              const double tmp77_1 = C_0_3*w35;
1934                              const register double tmp65_1 = C_1_2*w43;                              const double tmp65_1 = C_1_2*w43;
1935                              const register double tmp71_1 = C_0_3*w33;                              const double tmp71_1 = C_0_3*w33;
1936                              const register double tmp55_1 = C_1_1*w43;                              const double tmp55_1 = C_1_1*w43;
1937                              const register double tmp46_1 = tmp3_0*w28;                              const double tmp46_1 = tmp3_0*w28;
1938                              const register double tmp24_1 = C_0_0*w37;                              const double tmp24_1 = C_0_0*w37;
1939                              const register double tmp10_1 = tmp1_0*w39;                              const double tmp10_1 = tmp1_0*w39;
1940                              const register double tmp48_1 = C_1_0*w29;                              const double tmp48_1 = C_1_0*w29;
1941                              const register double tmp15_1 = C_0_1*w33;                              const double tmp15_1 = C_0_1*w33;
1942                              const register double tmp36_1 = C_1_2*w32;                              const double tmp36_1 = C_1_2*w32;
1943                              const register double tmp60_1 = C_1_0*w39;                              const double tmp60_1 = C_1_0*w39;
1944                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
1945                              const register double tmp16_1 = C_1_2*w29;                              const double tmp16_1 = C_1_2*w29;
1946                              const register double tmp1_1 = C_0_3*w37;                              const double tmp1_1 = C_0_3*w37;
1947                              const register double tmp7_1 = tmp2_0*w28;                              const double tmp7_1 = tmp2_0*w28;
1948                              const register double tmp39_1 = C_1_3*w40;                              const double tmp39_1 = C_1_3*w40;
1949                              const register double tmp44_1 = C_1_3*w42;                              const double tmp44_1 = C_1_3*w42;
1950                              const register double tmp57_1 = C_0_2*w38;                              const double tmp57_1 = C_0_2*w38;
1951                              const register double tmp78_1 = C_0_0*w34;                              const double tmp78_1 = C_0_0*w34;
1952                              const register double tmp11_1 = tmp0_0*w36;                              const double tmp11_1 = tmp0_0*w36;
1953                              const register double tmp23_1 = C_0_2*w33;                              const 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;  
1954                              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;
1955                              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;
1956                                EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1957                              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;
1958                              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;
1959                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1960                              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;
1961                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1962                              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;
1963                              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;  
1964                              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;
1965                              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;
1966                              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;
1967                              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;
1968                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1969                              const register double C_0 = C_p[0];                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1970                              const register double C_1 = C_p[1];                          } else { // constant data
1971                              const register double tmp1_1 = C_1*w45;                              const double C_0 = C_p[0];
1972                              const register double tmp3_1 = C_0*w51;                              const double C_1 = C_p[1];
1973                              const register double tmp4_1 = C_0*w44;                              const double tmp0_1 = C_0*w47;
1974                              const register double tmp7_1 = C_0*w46;                              const double tmp1_1 = C_1*w45;
1975                              const register double tmp5_1 = C_1*w49;                              const double tmp2_1 = C_1*w48;
1976                              const register double tmp2_1 = C_1*w48;                              const double tmp3_1 = C_0*w51;
1977                              const register double tmp0_1 = C_0*w47;                              const double tmp4_1 = C_0*w44;
1978                              const register double tmp6_1 = C_1*w50;                              const double tmp5_1 = C_1*w49;
1979                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const double tmp6_1 = C_1*w50;
1980                              EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;                              const double tmp7_1 = C_0*w46;
                             EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;  
1981                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1982                              EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
1983                                EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
1984                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
1985                              EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1986                                EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
1987                              EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;                              EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
1988                                EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
1989                              EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;                              EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
1990                              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;  
1991                              EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;                              EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
1992                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;                              EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
1993                              EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;                              EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
1994                              EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
1995                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
1996                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
1997                          }                          }
1998                      }                      }
1999                      ///////////////                      ///////////////
# Line 1883  void Rectangle::assemblePDESingle(Paso_S Line 2003  void Rectangle::assemblePDESingle(Paso_S
2003                          add_EM_S=true;                          add_EM_S=true;
2004                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2005                          if (D.actsExpanded()) {                          if (D.actsExpanded()) {
2006                              const register double D_0 = D_p[0];                              const double D_0 = D_p[0];
2007                              const register double D_1 = D_p[1];                              const double D_1 = D_p[1];
2008                              const register double D_2 = D_p[2];                              const double D_2 = D_p[2];
2009                              const register double D_3 = D_p[3];                              const double D_3 = D_p[3];
2010                              const register double tmp4_0 = D_1 + D_3;                              const double tmp0_0 = D_0 + D_1;
2011                              const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;                              const double tmp1_0 = D_2 + D_3;
2012                              const register double tmp5_0 = D_0 + D_2;                              const double tmp2_0 = D_0 + D_1 + D_2 + D_3;
2013                              const register double tmp0_0 = D_0 + D_1;                              const double tmp3_0 = D_1 + D_2;
2014                              const register double tmp6_0 = D_0 + D_3;                              const double tmp4_0 = D_1 + D_3;
2015                              const register double tmp1_0 = D_2 + D_3;                              const double tmp5_0 = D_0 + D_2;
2016                              const register double tmp3_0 = D_1 + D_2;                              const double tmp6_0 = D_0 + D_3;
2017                              const register double tmp16_1 = D_1*w56;                              const double tmp0_1 = tmp0_0*w52;
2018                              const register double tmp14_1 = tmp6_0*w54;                              const double tmp1_1 = tmp1_0*w53;
2019                              const register double tmp8_1 = D_3*w55;                              const double tmp2_1 = tmp2_0*w54;
2020                              const register double tmp2_1 = tmp2_0*w54;                              const double tmp3_1 = tmp1_0*w52;
2021                              const register double tmp12_1 = tmp5_0*w52;                              const double tmp4_1 = tmp0_0*w53;
2022                              const register double tmp4_1 = tmp0_0*w53;                              const double tmp5_1 = tmp3_0*w54;
2023                              const register double tmp3_1 = tmp1_0*w52;                              const double tmp6_1 = D_0*w55;
2024                              const register double tmp13_1 = tmp4_0*w53;                              const double tmp7_1 = D_3*w56;
2025                              const register double tmp10_1 = tmp4_0*w52;                              const double tmp8_1 = D_3*w55;
2026                              const register double tmp1_1 = tmp1_0*w53;                              const double tmp9_1 = D_0*w56;
2027                              const register double tmp7_1 = D_3*w56;                              const double tmp10_1 = tmp4_0*w52;
2028                              const register double tmp0_1 = tmp0_0*w52;                              const double tmp11_1 = tmp5_0*w53;
2029                              const register double tmp11_1 = tmp5_0*w53;                              const double tmp12_1 = tmp5_0*w52;
2030                              const register double tmp9_1 = D_0*w56;                              const double tmp13_1 = tmp4_0*w53;
2031                              const register double tmp5_1 = tmp3_0*w54;                              const double tmp14_1 = tmp6_0*w54;
2032                              const register double tmp18_1 = D_2*w56;                              const double tmp15_1 = D_2*w55;
2033                              const register double tmp17_1 = D_1*w55;                              const double tmp16_1 = D_1*w56;
2034                              const register double tmp6_1 = D_0*w55;                              const double tmp17_1 = D_1*w55;
2035                              const register double tmp15_1 = D_2*w55;                              const double tmp18_1 = D_2*w56;
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
2036                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2037                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2038                                EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
2039                              EM_S[INDEX2(3,0,4)]+=tmp2_1;                              EM_S[INDEX2(3,0,4)]+=tmp2_1;
2040                              EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2041                                EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
2042                              EM_S[INDEX2(2,1,4)]+=tmp2_1;                              EM_S[INDEX2(2,1,4)]+=tmp2_1;
2043                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
2044                              EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;                              EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
2045                              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;  
2046                              EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;                              EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
2047                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
2048                              EM_S[INDEX2(0,3,4)]+=tmp2_1;                              EM_S[INDEX2(0,3,4)]+=tmp2_1;
2049                              EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;                              EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
2050                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2051                              const register double D_0 = D_p[0];                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2052                              const register double tmp2_1 = D_0*w59;                          } else { // constant data
2053                              const register double tmp1_1 = D_0*w58;                              const double tmp0_1 = D_p[0]*w57;
2054                              const register double tmp0_1 = D_0*w57;                              const double tmp1_1 = D_p[0]*w58;
2055                              EM_S[INDEX2(0,1,4)]+=tmp0_1;                              const double tmp2_1 = D_p[0]*w59;
                             EM_S[INDEX2(1,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1;  
2056                              EM_S[INDEX2(0,0,4)]+=tmp2_1;                              EM_S[INDEX2(0,0,4)]+=tmp2_1;
2057                              EM_S[INDEX2(3,3,4)]+=tmp2_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
2058                                EM_S[INDEX2(2,0,4)]+=tmp0_1;
2059                              EM_S[INDEX2(3,0,4)]+=tmp1_1;                              EM_S[INDEX2(3,0,4)]+=tmp1_1;
2060                              EM_S[INDEX2(3,1,4)]+=tmp0_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1;
2061                                EM_S[INDEX2(1,1,4)]+=tmp2_1;
2062                              EM_S[INDEX2(2,1,4)]+=tmp1_1;                              EM_S[INDEX2(2,1,4)]+=tmp1_1;
2063                                EM_S[INDEX2(3,1,4)]+=tmp0_1;
2064                              EM_S[INDEX2(0,2,4)]+=tmp0_1;                              EM_S[INDEX2(0,2,4)]+=tmp0_1;
2065                              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;  
2066                              EM_S[INDEX2(2,2,4)]+=tmp2_1;                              EM_S[INDEX2(2,2,4)]+=tmp2_1;
2067                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(3,2,4)]+=tmp0_1;
2068                              EM_S[INDEX2(0,3,4)]+=tmp1_1;                              EM_S[INDEX2(0,3,4)]+=tmp1_1;
2069                              EM_S[INDEX2(1,1,4)]+=tmp2_1;                              EM_S[INDEX2(1,3,4)]+=tmp0_1;
2070                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
2071                                EM_S[INDEX2(3,3,4)]+=tmp2_1;
2072                          }                          }
2073                      }                      }
2074                      ///////////////                      ///////////////
# Line 1959  void Rectangle::assemblePDESingle(Paso_S Line 2078  void Rectangle::assemblePDESingle(Paso_S
2078                          add_EM_F=true;                          add_EM_F=true;
2079                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2080                          if (X.actsExpanded()) {                          if (X.actsExpanded()) {
2081                              const register double X_0_0 = X_p[INDEX2(0,0,2)];                              const double X_0_0 = X_p[INDEX2(0,0,2)];
2082                              const register double X_1_0 = X_p[INDEX2(1,0,2)];                              const double X_1_0 = X_p[INDEX2(1,0,2)];
2083                              const register double X_0_1 = X_p[INDEX2(0,1,2)];                              const double X_0_1 = X_p[INDEX2(0,1,2)];
2084                              const register double X_1_1 = X_p[INDEX2(1,1,2)];                              const double X_1_1 = X_p[INDEX2(1,1,2)];
2085                              const register double X_0_2 = X_p[INDEX2(0,2,2)];                              const double X_0_2 = X_p[INDEX2(0,2,2)];
2086                              const register double X_1_2 = X_p[INDEX2(1,2,2)];                              const double X_1_2 = X_p[INDEX2(1,2,2)];
2087                              const register double X_0_3 = X_p[INDEX2(0,3,2)];                              const double X_0_3 = X_p[INDEX2(0,3,2)];
2088                              const register double X_1_3 = X_p[INDEX2(1,3,2)];                              const double X_1_3 = X_p[INDEX2(1,3,2)];
2089                              const register double tmp1_0 = X_1_1 + X_1_3;                              const double tmp0_0 = X_0_2 + X_0_3;
2090                              const register double tmp3_0 = X_0_0 + X_0_1;                              const double tmp1_0 = X_1_1 + X_1_3;
2091                              const register double tmp2_0 = X_1_0 + X_1_2;                              const double tmp2_0 = X_1_0 + X_1_2;
2092                              const register double tmp0_0 = X_0_2 + X_0_3;                              const double tmp3_0 = X_0_0 + X_0_1;
2093                              const register double tmp8_1 = tmp2_0*w66;                              const double tmp0_1 = tmp0_0*w63;
2094                              const register double tmp5_1 = tmp3_0*w64;                              const double tmp1_1 = tmp1_0*w62;
2095                              const register double tmp14_1 = tmp0_0*w64;                              const double tmp2_1 = tmp2_0*w61;
2096                              const register double tmp3_1 = tmp3_0*w60;                              const double tmp3_1 = tmp3_0*w60;
2097                              const register double tmp9_1 = tmp3_0*w63;                              const double tmp4_1 = tmp0_0*w65;
2098                              const register double tmp13_1 = tmp3_0*w65;                              const double tmp5_1 = tmp3_0*w64;
2099                              const register double tmp12_1 = tmp1_0*w66;                              const double tmp6_1 = tmp2_0*w62;
2100                              const register double tmp10_1 = tmp0_0*w60;                              const double tmp7_1 = tmp1_0*w61;
2101                              const register double tmp2_1 = tmp2_0*w61;                              const double tmp8_1 = tmp2_0*w66;
2102                              const register double tmp6_1 = tmp2_0*w62;                              const double tmp9_1 = tmp3_0*w63;
2103                              const register double tmp4_1 = tmp0_0*w65;                              const double tmp10_1 = tmp0_0*w60;
2104                              const register double tmp11_1 = tmp1_0*w67;                              const double tmp11_1 = tmp1_0*w67;
2105                              const register double tmp1_1 = tmp1_0*w62;                              const double tmp12_1 = tmp1_0*w66;
2106                              const register double tmp7_1 = tmp1_0*w61;                              const double tmp13_1 = tmp3_0*w65;
2107                              const register double tmp0_1 = tmp0_0*w63;                              const double tmp14_1 = tmp0_0*w64;
2108                              const register double tmp15_1 = tmp2_0*w67;                              const double tmp15_1 = tmp2_0*w67;
2109                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2110                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2111                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2112                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2113                          } else { /* constant data */                          } else { // constant data
2114                              const register double X_0 = X_p[0];                              const double tmp0_1 = X_p[1]*w69;
2115                              const register double X_1 = X_p[1];                              const double tmp1_1 = X_p[0]*w68;
2116                              const register double tmp3_1 = X_1*w71;                              const double tmp2_1 = X_p[0]*w70;
2117                              const register double tmp2_1 = X_0*w70;                              const double tmp3_1 = X_p[1]*w71;
                             const register double tmp1_1 = X_0*w68;  
                             const register double tmp0_1 = X_1*w69;  
2118                              EM_F[0]+=tmp0_1 + tmp1_1;                              EM_F[0]+=tmp0_1 + tmp1_1;
2119                              EM_F[1]+=tmp0_1 + tmp2_1;                              EM_F[1]+=tmp0_1 + tmp2_1;
2120                              EM_F[2]+=tmp1_1 + tmp3_1;                              EM_F[2]+=tmp1_1 + tmp3_1;
# Line 2011  void Rectangle::assemblePDESingle(Paso_S Line 2128  void Rectangle::assemblePDESingle(Paso_S
2128                          add_EM_F=true;                          add_EM_F=true;
2129                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2130                          if (Y.actsExpanded()) {                          if (Y.actsExpanded()) {
2131                              const register double Y_0 = Y_p[0];                              const double Y_0 = Y_p[0];
2132                              const register double Y_1 = Y_p[1];                              const double Y_1 = Y_p[1];
2133                              const register double Y_2 = Y_p[2];                              const double Y_2 = Y_p[2];
2134                              const register double Y_3 = Y_p[3];                              const double Y_3 = Y_p[3];
2135                              const register double tmp0_0 = Y_1 + Y_2;                              const double tmp0_0 = Y_1 + Y_2;
2136                              const register double tmp1_0 = Y_0 + Y_3;                              const double tmp1_0 = Y_0 + Y_3;
2137                              const register double tmp9_1 = Y_0*w74;                              const double tmp0_1 = Y_0*w72;
2138                              const register double tmp4_1 = tmp1_0*w73;                              const double tmp1_1 = tmp0_0*w73;
2139                              const register double tmp5_1 = Y_2*w74;                              const double tmp2_1 = Y_3*w74;
2140                              const register double tmp7_1 = Y_1*w74;                              const double tmp3_1 = Y_1*w72;
2141                              const register double tmp6_1 = Y_2*w72;                              const double tmp4_1 = tmp1_0*w73;
2142                              const register double tmp2_1 = Y_3*w74;                              const double tmp5_1 = Y_2*w74;
2143                              const register double tmp8_1 = Y_3*w72;                              const double tmp6_1 = Y_2*w72;
2144                              const register double tmp3_1 = Y_1*w72;                              const double tmp7_1 = Y_1*w74;
2145                              const register double tmp0_1 = Y_0*w72;                              const double tmp8_1 = Y_3*w72;
2146                              const register double tmp1_1 = tmp0_0*w73;                              const double tmp9_1 = Y_0*w74;
2147                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2148                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2149                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2150                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2151                          } else { /* constant data */                          } else { // constant data
2152                              const register double Y_0 = Y_p[0];                              const double tmp0_1 = Y_p[0]*w75;
                             const register double tmp0_1 = Y_0*w75;  
2153                              EM_F[0]+=tmp0_1;                              EM_F[0]+=tmp0_1;
2154                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
2155                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
2156                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
2157                          }                          }
2158                      }                      }
                     /* GENERATOR SNIP_PDE_SINGLE BOTTOM */  
2159    
2160                      // 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)
2161                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
2162                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2163                      rowIndex.push_back(m_dofMap[firstNode]);                  } // end k0 loop
2164                      rowIndex.push_back(m_dofMap[firstNode+1]);              } // end k1 loop
2165                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);          } // end of colouring
2166                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);      } // end of parallel region
2167                      if (add_EM_F) {  }
2168                          //cout << "-- AddtoRHS -- " << endl;  
2169                          double *F_p=rhs.getSampleDataRW(0);  //protected
2170                          for (index_t i=0; i<rowIndex.size(); i++) {  void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2171                              if (rowIndex[i]<getNumDOF()) {          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2172                                  F_p[rowIndex[i]]+=EM_F[i];          const escript::Data& C, const escript::Data& D,
2173                                  //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;          const escript::Data& X, const escript::Data& Y) const
2174    {
2175        const double h0 = m_l0/m_gNE0;
2176        const double h1 = m_l1/m_gNE1;
2177        const double w0 = -.25*h1/h0;
2178        const double w1 = .25;
2179        const double w2 = -.25;
2180        const double w3 = .25*h0/h1;
2181        const double w4 = -.25*h0/h1;
2182        const double w5 = .25*h1/h0;
2183        const double w6 = -.125*h1;
2184        const double w7 = -.125*h0;
2185        const double w8 = .125*h1;
2186        const double w9 = .125*h0;
2187        const double w10 = .0625*h0*h1;
2188        const double w11 = -.5*h1;
2189        const double w12 = -.5*h0;
2190        const double w13 = .5*h1;
2191        const double w14 = .5*h0;
2192        const double w15 = .25*h0*h1;
2193    
2194        rhs.requireWrite();
2195    #pragma omp parallel
2196        {
2197            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2198    #pragma omp for
2199                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2200                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2201                        bool add_EM_S=false;
2202                        bool add_EM_F=false;
2203                        vector<double> EM_S(4*4, 0);
2204                        vector<double> EM_F(4, 0);
2205                        const index_t e = k0 + m_NE0*k1;
2206                        ///////////////
2207                        // process A //
2208                        ///////////////
2209                        if (!A.isEmpty()) {
2210                            add_EM_S=true;
2211                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2212                            const double A_00 = A_p[INDEX2(0,0,2)];
2213                            const double A_10 = A_p[INDEX2(1,0,2)];
2214                            const double A_01 = A_p[INDEX2(0,1,2)];
2215                            const double A_11 = A_p[INDEX2(1,1,2)];
2216                            const double tmp0_0 = A_01 + A_10;
2217                            const double tmp0_1 = A_11*w3;
2218                            const double tmp1_1 = A_00*w0;
2219                            const double tmp2_1 = A_01*w1;
2220                            const double tmp3_1 = A_10*w2;
2221                            const double tmp4_1 = tmp0_0*w1;
2222                            const double tmp5_1 = A_11*w4;
2223                            const double tmp6_1 = A_00*w5;
2224                            const double tmp7_1 = tmp0_0*w2;
2225                            const double tmp8_1 = A_10*w1;
2226                            const double tmp9_1 = A_01*w2;
2227                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2228                            EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2229                            EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2230                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2231                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2232                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2233                            EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2234                            EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2235                            EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2236                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2237                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2238                            EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2239                            EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2240                            EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2241                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2242                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2243                        }
2244                        ///////////////
2245                        // process B //
2246                        ///////////////
2247                        if (!B.isEmpty()) {
2248                            add_EM_S=true;
2249                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2250                            const double tmp2_1 = B_p[0]*w8;
2251                            const double tmp0_1 = B_p[0]*w6;
2252                            const double tmp3_1 = B_p[1]*w9;
2253                            const double tmp1_1 = B_p[1]*w7;
2254                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2255                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2256                            EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
2257                            EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;
2258                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2259                            EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;
2260                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;
2261                            EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;
2262                            EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;
2263                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
2264                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;
2265                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2266                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;
2267                            EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;
2268                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;
2269                            EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
2270                        }
2271                        ///////////////
2272                        // process C //
2273                        ///////////////
2274                        if (!C.isEmpty()) {
2275                            add_EM_S=true;
2276                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2277                            const double tmp1_1 = C_p[1]*w7;
2278                            const double tmp0_1 = C_p[0]*w8;
2279                            const double tmp3_1 = C_p[0]*w6;
2280                            const double tmp2_1 = C_p[1]*w9;
2281                            EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2282                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2283                            EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
2284                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
2285                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2286                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;
2287                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;
2288                            EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;
2289                            EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;
2290                            EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
2291                            EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
2292                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2293                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;
2294                            EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;
2295                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2296                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;
2297                        }
2298                        ///////////////
2299                        // process D //
2300                        ///////////////
2301                        if (!D.isEmpty()) {
2302                            add_EM_S=true;
2303                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2304                            const double tmp0_1 = D_p[0]*w10;
2305                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
2306                            EM_S[INDEX2(1,0,4)]+=tmp0_1;
2307                            EM_S[INDEX2(2,0,4)]+=tmp0_1;
2308                            EM_S[INDEX2(3,0,4)]+=tmp0_1;
2309                            EM_S[INDEX2(0,1,4)]+=tmp0_1;
2310                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
2311                            EM_S[INDEX2(2,1,4)]+=tmp0_1;
2312                            EM_S[INDEX2(3,1,4)]+=tmp0_1;
2313                            EM_S[INDEX2(0,2,4)]+=tmp0_1;
2314                            EM_S[INDEX2(1,2,4)]+=tmp0_1;
2315                            EM_S[INDEX2(2,2,4)]+=tmp0_1;
2316                            EM_S[INDEX2(3,2,4)]+=tmp0_1;
2317                            EM_S[INDEX2(0,3,4)]+=tmp0_1;
2318                            EM_S[INDEX2(1,3,4)]+=tmp0_1;
2319                            EM_S[INDEX2(2,3,4)]+=tmp0_1;
2320                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
2321                        }
2322                        ///////////////
2323                        // process X //
2324                        ///////////////
2325                        if (!X.isEmpty()) {
2326                            add_EM_F=true;
2327                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2328                            const double tmp0_1 = X_p[0]*w11;
2329                            const double tmp2_1 = X_p[0]*w13;
2330                            const double tmp1_1 = X_p[1]*w12;
2331                            const double tmp3_1 = X_p[1]*w14;
2332                            EM_F[0]+=tmp0_1 + tmp1_1;
2333                            EM_F[1]+=tmp1_1 + tmp2_1;
2334                            EM_F[2]+=tmp0_1 + tmp3_1;
2335                            EM_F[3]+=tmp2_1 + tmp3_1;
2336                        }
2337                        ///////////////
2338                        // process Y //
2339                        ///////////////
2340                        if (!Y.isEmpty()) {
2341                            add_EM_F=true;
2342                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2343                            const double tmp0_1 = Y_p[0]*w15;
2344                            EM_F[0]+=tmp0_1;
2345                            EM_F[1]+=tmp0_1;
2346                            EM_F[2]+=tmp0_1;
2347                            EM_F[3]+=tmp0_1;
2348                        }
2349    
2350                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2351                        const index_t firstNode=m_N0*k1+k0;
2352                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
2353                    } // end k0 loop
2354                } // end k1 loop
2355            } // end of colouring
2356        } // end of parallel region
2357    }
2358    
2359    //protected
2360    void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2361            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2362            const escript::Data& C, const escript::Data& D,
2363            const escript::Data& X, const escript::Data& Y) const
2364    {
2365        const double h0 = m_l0/m_gNE0;
2366        const double h1 = m_l1/m_gNE1;
2367        dim_t numEq, numComp;
2368        if (!mat)
2369            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
2370        else {
2371            numEq=mat->logical_row_block_size;
2372            numComp=mat->logical_col_block_size;
2373        }
2374    
2375        const double w0 = -0.1555021169820365539*h1/h0;
2376        const double w1 = 0.041666666666666666667;
2377        const double w2 = -0.15550211698203655390;
2378        const double w3 = 0.041666666666666666667*h0/h1;
2379        const double w4 = 0.15550211698203655390;
2380        const double w5 = -0.041666666666666666667;
2381        const double w6 = -0.01116454968463011277*h1/h0;
2382        const double w7 = 0.011164549684630112770;
2383        const double w8 = -0.011164549684630112770;
2384        const double w9 = -0.041666666666666666667*h1/h0;
2385        const double w10 = -0.041666666666666666667*h0/h1;
2386        const double w11 = 0.1555021169820365539*h1/h0;
2387        const double w12 = 0.1555021169820365539*h0/h1;
2388        const double w13 = 0.01116454968463011277*h0/h1;
2389        const double w14 = 0.01116454968463011277*h1/h0;
2390        const double w15 = 0.041666666666666666667*h1/h0;
2391        const double w16 = -0.01116454968463011277*h0/h1;
2392        const double w17 = -0.1555021169820365539*h0/h1;
2393        const double w18 = -0.33333333333333333333*h1/h0;
2394        const double w19 = 0.25000000000000000000;
2395        const double w20 = -0.25000000000000000000;
2396        const double w21 = 0.16666666666666666667*h0/h1;
2397        const double w22 = -0.16666666666666666667*h1/h0;
2398        const double w23 = -0.16666666666666666667*h0/h1;
2399        const double w24 = 0.33333333333333333333*h1/h0;
2400        const double w25 = 0.33333333333333333333*h0/h1;
2401        const double w26 = 0.16666666666666666667*h1/h0;
2402        const double w27 = -0.33333333333333333333*h0/h1;
2403        const double w28 = -0.032861463941450536761*h1;
2404        const double w29 = -0.032861463941450536761*h0;
2405        const double w30 = -0.12264065304058601714*h1;
2406        const double w31 = -0.0023593469594139828636*h1;
2407        const double w32 = -0.008805202725216129906*h0;
2408        const double w33 = -0.008805202725216129906*h1;
2409        const double w34 = 0.032861463941450536761*h1;
2410        const double w35 = 0.008805202725216129906*h1;
2411        const double w36 = 0.008805202725216129906*h0;
2412        const double w37 = 0.0023593469594139828636*h1;
2413        const double w38 = 0.12264065304058601714*h1;
2414        const double w39 = 0.032861463941450536761*h0;
2415        const double w40 = -0.12264065304058601714*h0;
2416        const double w41 = -0.0023593469594139828636*h0;
2417        const double w42 = 0.0023593469594139828636*h0;
2418        const double w43 = 0.12264065304058601714*h0;
2419        const double w44 = -0.16666666666666666667*h1;
2420        const double w45 = -0.083333333333333333333*h0;
2421        const double w46 = 0.083333333333333333333*h1;
2422        const double w47 = 0.16666666666666666667*h1;
2423        const double w48 = 0.083333333333333333333*h0;
2424        const double w49 = -0.16666666666666666667*h0;
2425        const double w50 = 0.16666666666666666667*h0;
2426        const double w51 = -0.083333333333333333333*h1;
2427        const double w52 = 0.025917019497006092316*h0*h1;
2428        const double w53 = 0.0018607582807716854616*h0*h1;
2429        const double w54 = 0.0069444444444444444444*h0*h1;
2430        const double w55 = 0.09672363354357992482*h0*h1;
2431        const double w56 = 0.00049858867864229740201*h0*h1;
2432        const double w57 = 0.055555555555555555556*h0*h1;
2433        const double w58 = 0.027777777777777777778*h0*h1;
2434        const double w59 = 0.11111111111111111111*h0*h1;
2435        const double w60 = -0.19716878364870322056*h1;
2436        const double w61 = -0.19716878364870322056*h0;
2437        const double w62 = -0.052831216351296779436*h0;
2438        const double w63 = -0.052831216351296779436*h1;
2439        const double w64 = 0.19716878364870322056*h1;
2440        const double w65 = 0.052831216351296779436*h1;
2441        const double w66 = 0.19716878364870322056*h0;
2442        const double w67 = 0.052831216351296779436*h0;
2443        const double w68 = -0.5*h1;
2444        const double w69 = -0.5*h0;
2445        const double w70 = 0.5*h1;
2446        const double w71 = 0.5*h0;
2447        const double w72 = 0.1555021169820365539*h0*h1;
2448        const double w73 = 0.041666666666666666667*h0*h1;
2449        const double w74 = 0.01116454968463011277*h0*h1;
2450        const double w75 = 0.25*h0*h1;
2451    
2452        rhs.requireWrite();
2453    #pragma omp parallel
2454        {
2455            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2456    #pragma omp for
2457                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2458                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2459                        bool add_EM_S=false;
2460                        bool add_EM_F=false;
2461                        vector<double> EM_S(4*4*numEq*numComp, 0);
2462                        vector<double> EM_F(4*numEq, 0);
2463                        const index_t e = k0 + m_NE0*k1;
2464                        ///////////////
2465                        // process A //
2466                        ///////////////
2467                        if (!A.isEmpty()) {
2468                            add_EM_S=true;
2469                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2470                            if (A.actsExpanded()) {
2471                                for (index_t k=0; k<numEq; k++) {
2472                                    for (index_t m=0; m<numComp; m++) {
2473                                        const double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];
2474                                        const double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];
2475                                        const double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];
2476                                        const double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];
2477                                        const double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];
2478                                        const double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];
2479                                        const double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];
2480                                        const double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];
2481                                        const double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];
2482                                        const double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];
2483                                        const double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];
2484                                        const double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];
2485                                        const double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];
2486                                        const double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];
2487                                        const double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];
2488                                        const double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];
2489                                        const double tmp0_0 = A_01_0 + A_01_3;
2490                                        const double tmp1_0 = A_00_0 + A_00_1;
2491                                        const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
2492                                        const double tmp3_0 = A_00_2 + A_00_3;
2493                                        const double tmp4_0 = A_10_1 + A_10_2;
2494                                        const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
2495                                        const double tmp6_0 = A_01_3 + A_10_0;
2496                                        const double tmp7_0 = A_01_0 + A_10_3;
2497                                        const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
2498                                        const double tmp9_0 = A_01_0 + A_10_0;
2499                                        const double tmp10_0 = A_01_3 + A_10_3;
2500                                        const double tmp11_0 = A_11_1 + A_11_3;
2501                                        const double tmp12_0 = A_11_0 + A_11_2;
2502                                        const double tmp13_0 = A_01_2 + A_10_1;
2503                                        const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
2504                                        const double tmp15_0 = A_01_1 + A_10_2;
2505                                        const double tmp16_0 = A_10_0 + A_10_3;
2506                                        const double tmp17_0 = A_01_1 + A_01_2;
2507                                        const double tmp18_0 = A_01_1 + A_10_1;
2508                                        const double tmp19_0 = A_01_2 + A_10_2;
2509                                        const double tmp0_1 = A_10_3*w8;
2510                                        const double tmp1_1 = tmp0_0*w1;
2511                                        const double tmp2_1 = A_01_1*w4;
2512                                        const double tmp3_1 = tmp1_0*w0;
2513                                        const double tmp4_1 = A_01_2*w7;
2514                                        const double tmp5_1 = tmp2_0*w3;
2515                                        const double tmp6_1 = tmp3_0*w6;
2516                                        const double tmp7_1 = A_10_0*w2;
2517                                        const double tmp8_1 = tmp4_0*w5;
2518                                        const double tmp9_1 = tmp2_0*w10;
2519                                        const double tmp10_1 = tmp5_0*w9;
2520                                        const double tmp11_1 = tmp6_0*w7;
2521                                        const double tmp12_1 = tmp7_0*w4;
2522                                        const double tmp13_1 = tmp8_0*w1;
2523                                        const double tmp14_1 = A_10_0*w8;
2524                                        const double tmp15_1 = A_01_2*w4;
2525                                        const double tmp16_1 = tmp3_0*w0;
2526                                        const double tmp17_1 = A_01_1*w7;
2527                                        const double tmp18_1 = tmp1_0*w6;
2528                                        const double tmp19_1 = A_10_3*w2;
2529                                        const double tmp20_1 = tmp9_0*w4;
2530                                        const double tmp21_1 = tmp1_0*w11;
2531                                        const double tmp22_1 = tmp10_0*w7;
2532                                        const double tmp23_1 = tmp3_0*w14;
2533                                        const double tmp24_1 = tmp11_0*w13;
2534                                        const double tmp25_1 = tmp12_0*w12;
2535                                        const double tmp26_1 = tmp10_0*w4;
2536                                        const double tmp27_1 = tmp3_0*w11;
2537                                        const double tmp28_1 = tmp9_0*w7;
2538                                        const double tmp29_1 = tmp1_0*w14;
2539                                        const double tmp30_1 = tmp12_0*w13;
2540                                        const double tmp31_1 = tmp11_0*w12;
2541                                        const double tmp32_1 = tmp13_0*w2;
2542                                        const double tmp33_1 = tmp14_0*w5;
2543                                        const double tmp34_1 = tmp15_0*w8;
2544                                        const double tmp35_1 = A_01_0*w8;
2545                                        const double tmp36_1 = tmp16_0*w1;
2546                                        const double tmp37_1 = A_10_1*w4;
2547                                        const double tmp38_1 = tmp5_0*w15;
2548                                        const double tmp39_1 = A_10_2*w7;
2549                                        const double tmp40_1 = tmp11_0*w17;
2550                                        const double tmp41_1 = A_01_3*w2;
2551                                        const double tmp42_1 = tmp12_0*w16;
2552                                        const double tmp43_1 = tmp17_0*w5;
2553                                        const double tmp44_1 = tmp7_0*w7;
2554                                        const double tmp45_1 = tmp6_0*w4;
2555                                        const double tmp46_1 = A_01_3*w8;
2556                                        const double tmp47_1 = A_10_2*w4;
2557                                        const double tmp48_1 = A_10_1*w7;
2558                                        const double tmp49_1 = tmp12_0*w17;
2559                                        const double tmp50_1 = A_01_0*w2;
2560                                        const double tmp51_1 = tmp11_0*w16;
2561                                        const double tmp52_1 = tmp18_0*w8;
2562                                        const double tmp53_1 = tmp19_0*w2;
2563                                        const double tmp54_1 = tmp13_0*w8;
2564                                        const double tmp55_1 = tmp15_0*w2;
2565                                        const double tmp56_1 = tmp19_0*w8;
2566                                        const double tmp57_1 = tmp18_0*w2;
2567                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
2568                                        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;
2569                                        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;
2570                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
2571                                        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;
2572                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
2573                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
2574                                        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;
2575                                        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;
2576                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
2577                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
2578                                        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;
2579                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
2580                                        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;
2581                                        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;
2582                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2583                                    }
2584                                }
2585                            } else { // constant data
2586                                for (index_t k=0; k<numEq; k++) {
2587                                    for (index_t m=0; m<numComp; m++) {
2588                                        const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
2589                                        const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
2590                                        const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
2591                                        const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
2592                                        const double tmp0_0 = A_01 + A_10;
2593                                        const double tmp0_1 = A_00*w18;
2594                                        const double tmp1_1 = A_01*w19;
2595                                        const double tmp2_1 = A_10*w20;
2596                                        const double tmp3_1 = A_11*w21;
2597                                        const double tmp4_1 = A_00*w22;
2598                                        const double tmp5_1 = tmp0_0*w19;
2599                                        const double tmp6_1 = A_11*w23;
2600                                        const double tmp7_1 = A_11*w25;
2601                                        const double tmp8_1 = A_00*w24;
2602                                        const double tmp9_1 = tmp0_0*w20;
2603                                        const double tmp10_1 = A_01*w20;
2604                                        const double tmp11_1 = A_11*w27;
2605                                        const double tmp12_1 = A_00*w26;
2606                                        const double tmp13_1 = A_10*w19;
2607                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2608                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2609                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
2610                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2611                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2612                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2613                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
2614                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
2615                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
2616                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
2617                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2618                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2619                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2620                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
2621                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2622                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2623                                    }
2624                              }                              }
2625                          }                          }
                         //cout << "---"<<endl;  
2626                      }                      }
2627                      if (add_EM_S) {                      ///////////////
2628                          //cout << "-- AddtoSystem -- " << endl;                      // process B //
2629                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);                      ///////////////
2630                        if (!B.isEmpty()) {
2631                            add_EM_S=true;
2632                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2633                            if (B.actsExpanded()) {
2634                                for (index_t k=0; k<numEq; k++) {
2635                                    for (index_t m=0; m<numComp; m++) {
2636                                        const double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
2637                                        const double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
2638                                        const double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
2639                                        const double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
2640                                        const double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];
2641                                        const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
2642                                        const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
2643                                        const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
2644                                        const double tmp0_0 = B_1_0 + B_1_1;
2645                                        const double tmp1_0 = B_1_2 + B_1_3;
2646                                        const double tmp2_0 = B_0_1 + B_0_3;
2647                                        const double tmp3_0 = B_0_0 + B_0_2;
2648                                        const double tmp63_1 = B_1_1*w42;
2649                                        const double tmp79_1 = B_1_1*w40;
2650                                        const double tmp37_1 = tmp3_0*w35;
2651                                        const double tmp8_1 = tmp0_0*w32;
2652                                        const double tmp71_1 = B_0_1*w34;
2653                                        const double tmp19_1 = B_0_3*w31;
2654                                        const double tmp15_1 = B_0_3*w34;
2655                                        const double tmp9_1 = tmp3_0*w34;
2656                                        const double tmp35_1 = B_1_0*w36;
2657                                        const double tmp66_1 = B_0_3*w28;
2658                                        const double tmp28_1 = B_1_0*w42;
2659                                        const double tmp22_1 = B_1_0*w40;
2660                                        const double tmp16_1 = B_1_2*w29;
2661                                        const double tmp6_1 = tmp2_0*w35;
2662                                        const double tmp55_1 = B_1_3*w40;
2663                                        const double tmp50_1 = B_1_3*w42;
2664                                        const double tmp7_1 = tmp1_0*w29;
2665                                        const double tmp1_1 = tmp1_0*w32;
2666                                        const double tmp57_1 = B_0_3*w30;
2667                                        const double tmp18_1 = B_1_1*w32;
2668                                        const double tmp53_1 = B_1_0*w41;
2669                                        const double tmp61_1 = B_1_3*w36;
2670                                        const double tmp27_1 = B_0_3*w38;
2671                                        const double tmp64_1 = B_0_2*w30;
2672                                        const double tmp76_1 = B_0_1*w38;
2673                                        const double tmp39_1 = tmp2_0*w34;
2674                                        const double tmp62_1 = B_0_1*w31;
2675                                        const double tmp56_1 = B_0_0*w31;
2676                                        const double tmp49_1 = B_1_1*w36;
2677                                        const double tmp2_1 = B_0_2*w31;
2678                                        const double tmp23_1 = B_0_2*w33;
2679                                        const double tmp38_1 = B_1_1*w43;
2680                                        const double tmp74_1 = B_1_2*w41;
2681                                        const double tmp43_1 = B_1_1*w41;
2682                                        const double tmp58_1 = B_0_2*w28;
2683                                        const double tmp67_1 = B_0_0*w33;
2684                                        const double tmp33_1 = tmp0_0*w39;
2685                                        const double tmp4_1 = B_0_0*w28;
2686                                        const double tmp20_1 = B_0_0*w30;
2687                                        const double tmp13_1 = B_0_2*w38;
2688                                        const double tmp65_1 = B_1_2*w43;
2689                                        const double tmp0_1 = tmp0_0*w29;
2690                                        const double tmp41_1 = tmp3_0*w33;
2691                                        const double tmp73_1 = B_0_2*w37;
2692                                        const double tmp69_1 = B_0_0*w38;
2693                                        const double tmp48_1 = B_1_2*w39;
2694                                        const double tmp59_1 = B_0_1*w33;
2695                                        const double tmp17_1 = B_1_3*w41;
2696                                        const double tmp5_1 = B_0_3*w33;
2697                                        const double tmp3_1 = B_0_1*w30;
2698                                        const double tmp21_1 = B_0_1*w28;
2699                                        const double tmp42_1 = B_1_0*w29;
2700                                        const double tmp54_1 = B_1_2*w32;
2701                                        const double tmp60_1 = B_1_0*w39;
2702                                        const double tmp32_1 = tmp1_0*w36;
2703                                        const double tmp10_1 = B_0_1*w37;
2704                                        const double tmp14_1 = B_0_0*w35;
2705                                        const double tmp29_1 = B_0_1*w35;
2706                                        const double tmp26_1 = B_1_2*w36;
2707                                        const double tmp30_1 = B_1_3*w43;
2708                                        const double tmp70_1 = B_0_2*w35;
2709                                        const double tmp34_1 = B_1_3*w39;
2710                                        const double tmp51_1 = B_1_0*w43;
2711                                        const double tmp31_1 = B_0_2*w34;
2712                                        const double tmp45_1 = tmp3_0*w28;
2713                                        const double tmp11_1 = tmp1_0*w39;
2714                                        const double tmp52_1 = B_1_1*w29;
2715                                        const double tmp44_1 = B_1_3*w32;
2716                                        const double tmp25_1 = B_1_1*w39;
2717                                        const double tmp47_1 = tmp2_0*w33;
2718                                        const double tmp72_1 = B_1_3*w29;
2719                                        const double tmp40_1 = tmp2_0*w28;
2720                                        const double tmp46_1 = B_1_2*w40;
2721                                        const double tmp36_1 = B_1_2*w42;
2722                                        const double tmp24_1 = B_0_0*w37;
2723                                        const double tmp77_1 = B_0_3*w35;
2724                                        const double tmp68_1 = B_0_3*w37;
2725                                        const double tmp78_1 = B_0_0*w34;
2726                                        const double tmp12_1 = tmp0_0*w36;
2727                                        const double tmp75_1 = B_1_0*w32;
2728                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2729                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2730                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2731                                        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;
2732                                        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;
2733                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
2734                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
2735                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
2736                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
2737                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2738                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
2739                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2740                                        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;
2741                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2742                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
2743                                        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;
2744                                    }
2745                                }
2746                            } else { // constant data
2747                                for (index_t k=0; k<numEq; k++) {
2748                                    for (index_t m=0; m<numComp; m++) {
2749                                        const double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
2750                                        const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
2751                                        const double tmp6_1 = B_1*w50;
2752                                        const double tmp1_1 = B_1*w45;
2753                                        const double tmp5_1 = B_1*w49;
2754                                        const double tmp4_1 = B_1*w48;
2755                                        const double tmp0_1 = B_0*w44;
2756                                        const double tmp2_1 = B_0*w46;
2757                                        const double tmp7_1 = B_0*w51;
2758                                        const double tmp3_1 = B_0*w47;
2759                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
2760                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
2761                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
2762                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;
2763                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;
2764                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
2765                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;
2766                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;
2767                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;
2768                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;
2769                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;
2770                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;
2771                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;
2772                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
2773                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;
2774                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;
2775                                    }
2776                                }
2777                            }
2778                        }
2779                        ///////////////
2780                        // process C //
2781                        ///////////////
2782                        if (!C.isEmpty()) {
2783                            add_EM_S=true;
2784                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2785                            if (C.actsExpanded()) {
2786                                for (index_t k=0; k<numEq; k++) {
2787                                    for (index_t m=0; m<numComp; m++) {
2788                                        const double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
2789                                        const double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
2790                                        const double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
2791                                        const double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
2792                                        const double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];
2793                                        const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
2794                                        const double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];
2795                                        const double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];
2796                                        const double tmp2_0 = C_0_1 + C_0_3;
2797                                        const double tmp1_0 = C_1_2 + C_1_3;
2798                                        const double tmp3_0 = C_0_0 + C_0_2;
2799                                        const double tmp0_0 = C_1_0 + C_1_1;
2800                                        const double tmp64_1 = C_0_2*w30;
2801                                        const double tmp14_1 = C_0_2*w28;
2802                                        const double tmp19_1 = C_0_3*w31;
2803                                        const double tmp22_1 = C_1_0*w40;
2804                                        const double tmp37_1 = tmp3_0*w35;
2805                                        const double tmp29_1 = C_0_1*w35;
2806                                        const double tmp73_1 = C_0_2*w37;
2807                                        const double tmp74_1 = C_1_2*w41;
2808                                        const double tmp52_1 = C_1_3*w39;
2809                                        const double tmp25_1 = C_1_1*w39;
2810                                        const double tmp62_1 = C_0_1*w31;
2811                                        const double tmp79_1 = C_1_1*w40;
2812                                        const double tmp43_1 = C_1_1*w36;
2813                                        const double tmp27_1 = C_0_3*w38;
2814                                        const double tmp28_1 = C_1_0*w42;
2815                                        const double tmp63_1 = C_1_1*w42;
2816                                        const double tmp59_1 = C_0_3*w34;
2817                                        const double tmp72_1 = C_1_3*w29;
2818                                        const double tmp40_1 = tmp2_0*w35;
2819                                        const double tmp13_1 = C_0_3*w30;
2820                                        const double tmp51_1 = C_1_2*w40;
2821                                        const double tmp54_1 = C_1_2*w42;
2822                                        const double tmp12_1 = C_0_0*w31;
2823                                        const double tmp2_1 = tmp1_0*w32;
2824                                        const double tmp68_1 = C_0_2*w31;
2825                                        const double tmp75_1 = C_1_0*w32;
2826                                        const double tmp49_1 = C_1_1*w41;
2827                                        const double tmp4_1 = C_0_2*w35;
2828                                        const double tmp66_1 = C_0_3*w28;
2829                                        const double tmp56_1 = C_0_1*w37;
2830                                        const double tmp5_1 = C_0_1*w34;
2831                                        const double tmp38_1 = tmp2_0*w34;
2832                                        const double tmp76_1 = C_0_1*w38;
2833                                        const double tmp21_1 = C_0_1*w28;
2834                                        const double tmp69_1 = C_0_1*w30;
2835                                        const double tmp53_1 = C_1_0*w36;
2836                                        const double tmp42_1 = C_1_2*w39;
2837                                        const double tmp32_1 = tmp1_0*w29;
2838                                        const double tmp45_1 = C_1_0*w43;
2839                                        const double tmp33_1 = tmp0_0*w32;
2840                                        const double tmp35_1 = C_1_0*w41;
2841                                        const double tmp26_1 = C_1_2*w36;
2842                                        const double tmp67_1 = C_0_0*w33;
2843                                        const double tmp31_1 = C_0_2*w34;
2844                                        const double tmp20_1 = C_0_0*w30;
2845                                        const double tmp70_1 = C_0_0*w28;
2846                                        const double tmp8_1 = tmp0_0*w39;
2847                                        const double tmp30_1 = C_1_3*w43;
2848