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

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

Parent Directory | Revision Log | View Patch Patch

-revision 3760 by caltinay,
Mon Jan  9 05:21:18 2012 UTC
+revision 3777 by caltinay,
Thu Jan 19 06:17:38 2012 UTC
 Line 1
  /*******************************************************
  *
- * Copyright (c) 2003-2011 by University of Queensland
+ * Copyright (c) 2003-2012 by University of Queensland
  * Earth Systems Science Computational Center (ESSCC)
  * http://www.uq.edu.au/esscc
  *
 Line 49 
 Rectangle::Rectangle(int n0, int n1, dou
      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
          throw RipleyException("Too few elements for the number of ranks");
-     // local number of elements (including overlap)
+     // local number of elements (with and without overlap)
-     m_NE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
+     m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
          m_NE0++;
-     m_NE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
+     else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)
+         m_ownNE0--;
+     m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)
          m_NE1++;
+     else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)
+         m_ownNE1--;
      // local number of nodes
      m_N0 = m_NE0+1;
-Line 106 
 void Rectangle::dump(const string& fileN
+Line 111 
 void Rectangle::dump(const string& fileN
      const int NUM_SILO_FILES = 1;
      const char* blockDirFmt = "/block%04d";
      int driver=DB_HDF5;
-     string siloPath;
      DBfile* dbfile = NULL;
  #ifdef ESYS_MPI
-Line 126 
 void Rectangle::dump(const string& fileN
+Line 130 
 void Rectangle::dump(const string& fileN
                          PMPIO_DefaultClose, (void*)&driver);
          }
          if (baton) {
-             char str[64];
+             char siloPath[64];
-             snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
+             snprintf(siloPath, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
-             siloPath = str;
+             dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath);
-             dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
          }
  #endif
      } else {
-Line 141 
 void Rectangle::dump(const string& fileN
+Line 144 
 void Rectangle::dump(const string& fileN
              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
                      getDescription().c_str(), driver);
          }
+         char siloPath[64];
+         snprintf(siloPath, 64, blockDirFmt, 0);
+         DBMkDir(dbfile, siloPath);
+         DBSetDir(dbfile, siloPath);
      }
      if (!dbfile)
-Line 241 
 const int* Rectangle::borrowSampleRefere
+Line 248 
 const int* Rectangle::borrowSampleRefere
  {
      switch (fsType) {
          case Nodes:
-         case ReducedNodes: //FIXME: reduced
+         case ReducedNodes: // FIXME: reduced
              return &m_nodeId[0];
          case DegreesOfFreedom:
-         case ReducedDegreesOfFreedom: //FIXME: reduced
+         case ReducedDegreesOfFreedom: // FIXME: reduced
              return &m_dofId[0];
          case Elements:
          case ReducedElements:
-Line 269 
 bool Rectangle::ownSample(int fsType, in
+Line 276 
 bool Rectangle::ownSample(int fsType, in
      switch (fsType) {
          case Nodes:
-         case ReducedNodes: //FIXME: reduced
+         case ReducedNodes: // FIXME: reduced
              return (m_dofMap[id] < getNumDOF());
          case DegreesOfFreedom:
          case ReducedDegreesOfFreedom:
-Line 281 
 bool Rectangle::ownSample(int fsType, in
+Line 288 
 bool Rectangle::ownSample(int fsType, in
          case FaceElements:
          case ReducedFaceElements:
              {
-                 // check ownership of face element's first node
+                 // determine which face the sample belongs to before
+                 // checking ownership of corresponding element's first node
                  const IndexVector faces = getNumFacesPerBoundary();
                  dim_t n=0;
                  for (size_t i=0; i<faces.size(); i++) {
-Line 289 
 bool Rectangle::ownSample(int fsType, in
+Line 297 
 bool Rectangle::ownSample(int fsType, in
                      if (id<n) {
                          index_t k;
                          if (i==1)
-                             k=m_N0-1;
+                             k=m_N0-2;
                          else if (i==3)
-                             k=m_N0*(m_N1-1);
+                             k=m_N0*(m_N1-2);
                          else
                              k=0;
                          // determine whether to move right or up
-Line 311 
 bool Rectangle::ownSample(int fsType, in
+Line 319 
 bool Rectangle::ownSample(int fsType, in
      throw RipleyException(msg.str());
  }
- void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const
- {
-     escript::Data& in = *const_cast<escript::Data*>(&cIn);
-     const dim_t numComp = in.getDataPointSize();
-     const double h0 = m_l0/m_gNE0;
-     const double h1 = m_l1/m_gNE1;
-     const double cx0 = -1./h0;
-     const double cx1 = -.78867513459481288225/h0;
-     const double cx2 = -.5/h0;
-     const double cx3 = -.21132486540518711775/h0;
-     const double cx4 = .21132486540518711775/h0;
-     const double cx5 = .5/h0;
-     const double cx6 = .78867513459481288225/h0;
-     const double cx7 = 1./h0;
-     const double cy0 = -1./h1;
-     const double cy1 = -.78867513459481288225/h1;
-     const double cy2 = -.5/h1;
-     const double cy3 = -.21132486540518711775/h1;
-     const double cy4 = .21132486540518711775/h1;
-     const double cy5 = .5/h1;
-     const double cy6 = .78867513459481288225/h1;
-     const double cy7 = 1./h1;
-     if (out.getFunctionSpace().getTypeCode() == Elements) {
-         out.requireWrite();
-         /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */
- #pragma omp parallel for
-         for (index_t k1=0; k1 < m_NE1; ++k1) {
-             for (index_t k0=0; k0 < m_NE0; ++k0) {
-                 const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
-                 const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
-                 const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
-                 const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
-                 double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
-                 for (index_t i=0; i < numComp; ++i) {
-                     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,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
-                     o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
-                     o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
-                     o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
-                     o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
-                     o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
-                     o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
-                 } /* end of component loop i */
-             } /* end of k0 loop */
-         } /* end of k1 loop */
-         /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
-     } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
-         out.requireWrite();
-         /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */
- #pragma omp parallel for
-         for (index_t k1=0; k1 < m_NE1; ++k1) {
-             for (index_t k0=0; k0 < m_NE0; ++k0) {
-                 const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
-                 const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
-                 const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
-                 const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
-                 double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
-                 for (index_t i=0; i < numComp; ++i) {
-                     o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
-                     o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
-                 } /* end of component loop i */
-             } /* end of k0 loop */
-         } /* end of k1 loop */
-         /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
-     } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
-         out.requireWrite();
- #pragma omp parallel
-         {
-             /*** GENERATOR SNIP_GRAD_FACES TOP */
-             if (m_faceOffset[0] > -1) {
- #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
-                     const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
-                     const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
-                     const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
-                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
-                     for (index_t i=0; i < numComp; ++i) {
-                         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,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
-                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
-                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
-                     } /* end of component loop i */
-                 } /* end of k1 loop */
-             } /* end of face 0 */
-             if (m_faceOffset[1] > -1) {
- #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
-                     const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
-                     const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
-                     const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
-                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
-                     for (index_t i=0; i < numComp; ++i) {
-                         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,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
-                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
-                         o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
-                     } /* end of component loop i */
-                 } /* end of k1 loop */
-             } /* end of face 1 */
-             if (m_faceOffset[2] > -1) {
- #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
-                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
-                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
-                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
-                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
-                     for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
-                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
-                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
-                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
-                     } /* end of component loop i */
-                 } /* end of k0 loop */
-             } /* end of face 2 */
-             if (m_faceOffset[3] > -1) {
- #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
-                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
-                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
-                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
-                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
-                     for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
-                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
-                         o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
-                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
-                     } /* end of component loop i */
-                 } /* end of k0 loop */
-             } /* end of face 3 */
-             /* GENERATOR SNIP_GRAD_FACES BOTTOM */
-         } // end of parallel section
-     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
-         out.requireWrite();
- #pragma omp parallel
-         {
-             /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
-             if (m_faceOffset[0] > -1) {
- #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
-                     const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
-                     const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
-                     const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
-                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
-                     for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
-                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
-                     } /* end of component loop i */
-                 } /* end of k1 loop */
-             } /* end of face 0 */
-             if (m_faceOffset[1] > -1) {
- #pragma omp for nowait
-                 for (index_t k1=0; k1 < m_NE1; ++k1) {
-                     const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
-                     const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
-                     const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
-                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
-                     for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
-                         o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
-                     } /* end of component loop i */
-                 } /* end of k1 loop */
-             } /* end of face 1 */
-             if (m_faceOffset[2] > -1) {
- #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
-                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
-                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
-                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
-                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
-                     for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
-                         o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
-                     } /* end of component loop i */
-                 } /* end of k0 loop */
-             } /* end of face 2 */
-             if (m_faceOffset[3] > -1) {
- #pragma omp for nowait
-                 for (index_t k0=0; k0 < m_NE0; ++k0) {
-                     const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
-                     const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
-                     const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
-                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
-                     for (index_t i=0; i < numComp; ++i) {
-                         o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
-                         o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
-                     } /* end of component loop i */
-                 } /* end of k0 loop */
-             } /* end of face 3 */
-             /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */
-         } // end of parallel section
-     } else {
-         stringstream msg;
-         msg << "setToGradient() not implemented for "
-             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
-         throw RipleyException(msg.str());
-     }
- }
- void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
- {
-     escript::Data& in = *const_cast<escript::Data*>(&arg);
-     const dim_t numComp = in.getDataPointSize();
-     const double h0 = m_l0/m_gNE0;
-     const double h1 = m_l1/m_gNE1;
-     if (arg.getFunctionSpace().getTypeCode() == Elements) {
-         const double w_0 = h0*h1/4.;
- #pragma omp parallel
-         {
-             vector<double> int_local(numComp, 0);
- #pragma omp for nowait
-             for (index_t k1 = 0; k1 < m_NE1; ++k1) {
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
-                     const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
-                     for (index_t i=0; i < numComp; ++i) {
-                         const register double f_0 = f[INDEX2(i,0,numComp)];
-                         const register double f_1 = f[INDEX2(i,1,numComp)];
-                         const register double f_2 = f[INDEX2(i,2,numComp)];
-                         const register double f_3 = f[INDEX2(i,3,numComp)];
-                         int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
-                     }  /* end of component loop i */
-                 } /* end of k0 loop */
-             } /* end of k1 loop */
- #pragma omp critical
-             for (index_t i=0; i<numComp; i++)
-                 integrals[i]+=int_local[i];
-         } // end of parallel section
-     } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
-         const double w_0 = h0*h1;
- #pragma omp parallel
-         {
-             vector<double> int_local(numComp, 0);
- #pragma omp for nowait
-             for (index_t k1 = 0; k1 < m_NE1; ++k1) {
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
-                     const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
-                     for (index_t i=0; i < numComp; ++i) {
-                         int_local[i]+=f[i]*w_0;
-                     }  /* end of component loop i */
-                 } /* end of k0 loop */
-             } /* end of k1 loop */
- #pragma omp critical
-             for (index_t i=0; i<numComp; i++)
-                 integrals[i]+=int_local[i];
-         } // end of parallel section
-     } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
-         const double w_0 = h0/2.;
-         const double w_1 = h1/2.;
- #pragma omp parallel
-         {
-             vector<double> int_local(numComp, 0);
-             if (m_faceOffset[0] > -1) {
- #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
-                     const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
-                     for (index_t i=0; i < numComp; ++i) {
-                         const register double f_0 = f[INDEX2(i,0,numComp)];
-                         const register double f_1 = f[INDEX2(i,1,numComp)];
-                         int_local[i]+=(f_0+f_1)*w_1;
-                     }  /* end of component loop i */
-                 } /* end of k1 loop */
-             }
-             if (m_faceOffset[1] > -1) {
- #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
-                     const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
-                     for (index_t i=0; i < numComp; ++i) {
-                         const register double f_0 = f[INDEX2(i,0,numComp)];
-                         const register double f_1 = f[INDEX2(i,1,numComp)];
-                         int_local[i]+=(f_0+f_1)*w_1;
-                     }  /* end of component loop i */
-                 } /* end of k1 loop */
-             }
-             if (m_faceOffset[2] > -1) {
- #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
-                     const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
-                     for (index_t i=0; i < numComp; ++i) {
-                         const register double f_0 = f[INDEX2(i,0,numComp)];
-                         const register double f_1 = f[INDEX2(i,1,numComp)];
-                         int_local[i]+=(f_0+f_1)*w_0;
-                     }  /* end of component loop i */
-                 } /* end of k0 loop */
-             }
-             if (m_faceOffset[3] > -1) {
- #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
-                     const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
-                     for (index_t i=0; i < numComp; ++i) {
-                         const register double f_0 = f[INDEX2(i,0,numComp)];
-                         const register double f_1 = f[INDEX2(i,1,numComp)];
-                         int_local[i]+=(f_0+f_1)*w_0;
-                     }  /* end of component loop i */
-                 } /* end of k0 loop */
-             }
- #pragma omp critical
-             for (index_t i=0; i<numComp; i++)
-                 integrals[i]+=int_local[i];
-         } // end of parallel section
-     } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
- #pragma omp parallel
-         {
-             vector<double> int_local(numComp, 0);
-             if (m_faceOffset[0] > -1) {
- #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
-                     const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
-                     for (index_t i=0; i < numComp; ++i) {
-                         int_local[i]+=f[i]*h1;
-                     }  /* end of component loop i */
-                 } /* end of k1 loop */
-             }
-             if (m_faceOffset[1] > -1) {
- #pragma omp for nowait
-                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {
-                     const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
-                     for (index_t i=0; i < numComp; ++i) {
-                         int_local[i]+=f[i]*h1;
-                     }  /* end of component loop i */
-                 } /* end of k1 loop */
-             }
-             if (m_faceOffset[2] > -1) {
- #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
-                     const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
-                     for (index_t i=0; i < numComp; ++i) {
-                         int_local[i]+=f[i]*h0;
-                     }  /* end of component loop i */
-                 } /* end of k0 loop */
-             }
-             if (m_faceOffset[3] > -1) {
- #pragma omp for nowait
-                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {
-                     const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
-                     for (index_t i=0; i < numComp; ++i) {
-                         int_local[i]+=f[i]*h0;
-                     }  /* end of component loop i */
-                 } /* end of k0 loop */
-             }
- #pragma omp critical
-             for (index_t i=0; i<numComp; i++)
-                 integrals[i]+=int_local[i];
-         } // 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) {
-Line 849 
 void Rectangle::setToSize(escript::Data&
+Line 488 
 void Rectangle::setToSize(escript::Data&
  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
                                                  bool reducedColOrder) const
  {
+     /* FIXME: reduced
      if (reducedRowOrder || reducedColOrder)
          throw RipleyException("getPattern() not implemented for reduced order");
+     */
      return m_pattern;
  }
-Line 906 
 IndexVector Rectangle::getNumFacesPerBou
+Line 546 
 IndexVector Rectangle::getNumFacesPerBou
      return ret;
  }
+ IndexVector Rectangle::getNumSubdivisionsPerDim() const
+ {
+     IndexVector ret;
+     ret.push_back(m_NX);
+     ret.push_back(m_NY);
+     return ret;
+ }
  pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
  {
      if (dim==0) {
-Line 956 
 void Rectangle::assembleCoordinates(escr
+Line 604 
 void Rectangle::assembleCoordinates(escr
  }
  //protected
+ void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
+ {
+     const dim_t numComp = in.getDataPointSize();
+     const double h0 = m_l0/m_gNE0;
+     const double h1 = m_l1/m_gNE1;
+     const double cx0 = -1./h0;
+     const double cx1 = -.78867513459481288225/h0;
+     const double cx2 = -.5/h0;
+     const double cx3 = -.21132486540518711775/h0;
+     const double cx4 = .21132486540518711775/h0;
+     const double cx5 = .5/h0;
+     const double cx6 = .78867513459481288225/h0;
+     const double cx7 = 1./h0;
+     const double cy0 = -1./h1;
+     const double cy1 = -.78867513459481288225/h1;
+     const double cy2 = -.5/h1;
+     const double cy3 = -.21132486540518711775/h1;
+     const double cy4 = .21132486540518711775/h1;
+     const double cy5 = .5/h1;
+     const double cy6 = .78867513459481288225/h1;
+     const double cy7 = 1./h1;
+     if (out.getFunctionSpace().getTypeCode() == Elements) {
+         out.requireWrite();
+ #pragma omp parallel for
+         for (index_t k1=0; k1 < m_NE1; ++k1) {
+             for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
+                 const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
+                 const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
+                 const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
+                 double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
+                 for (index_t i=0; i < numComp; ++i) {
+                     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,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
+                     o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
+                     o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
+                     o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
+                     o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
+                     o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
+                     o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
+                 } /* end of component loop i */
+             } /* end of k0 loop */
+         } /* end of k1 loop */
+     } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
+         out.requireWrite();
+ #pragma omp parallel for
+         for (index_t k1=0; k1 < m_NE1; ++k1) {
+             for (index_t k0=0; k0 < m_NE0; ++k0) {
+                 const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
+                 const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
+                 const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
+                 const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
+                 double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
+                 for (index_t i=0; i < numComp; ++i) {
+                     o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
+                     o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
+                 } /* end of component loop i */
+             } /* end of k0 loop */
+         } /* end of k1 loop */
+     } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
+         out.requireWrite();
+ #pragma omp parallel
+         {
+             if (m_faceOffset[0] > -1) {
+ #pragma omp for nowait
+                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                     const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
+                     const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
+                     const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
+                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
+                     for (index_t i=0; i < numComp; ++i) {
+                         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,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
+                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
+                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
+                     } /* end of component loop i */
+                 } /* end of k1 loop */
+             } /* end of face 0 */
+             if (m_faceOffset[1] > -1) {
+ #pragma omp for nowait
+                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                     const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
+                     const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
+                     const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
+                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
+                     for (index_t i=0; i < numComp; ++i) {
+                         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,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
+                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
+                         o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
+                     } /* end of component loop i */
+                 } /* end of k1 loop */
+             } /* end of face 1 */
+             if (m_faceOffset[2] > -1) {
+ #pragma omp for nowait
+                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                     const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
+                     const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
+                     const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
+                     const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
+                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
+                     for (index_t i=0; i < numComp; ++i) {
+                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
+                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
+                         o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
+                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
+                     } /* end of component loop i */
+                 } /* end of k0 loop */
+             } /* end of face 2 */
+             if (m_faceOffset[3] > -1) {
+ #pragma omp for nowait
+                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                     const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
+                     const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
+                     const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
+                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
+                     for (index_t i=0; i < numComp; ++i) {
+                         o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
+                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
+                         o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
+                         o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
+                     } /* end of component loop i */
+                 } /* end of k0 loop */
+             } /* end of face 3 */
+         } // end of parallel section
+     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
+         out.requireWrite();
+ #pragma omp parallel
+         {
+             if (m_faceOffset[0] > -1) {
+ #pragma omp for nowait
+                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                     const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
+                     const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
+                     const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
+                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
+                     for (index_t i=0; i < numComp; ++i) {
+                         o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
+                         o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
+                     } /* end of component loop i */
+                 } /* end of k1 loop */
+             } /* end of face 0 */
+             if (m_faceOffset[1] > -1) {
+ #pragma omp for nowait
+                 for (index_t k1=0; k1 < m_NE1; ++k1) {
+                     const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
+                     const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
+                     const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
+                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
+                     for (index_t i=0; i < numComp; ++i) {
+                         o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
+                         o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
+                     } /* end of component loop i */
+                 } /* end of k1 loop */
+             } /* end of face 1 */
+             if (m_faceOffset[2] > -1) {
+ #pragma omp for nowait
+                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                     const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
+                     const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
+                     const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
+                     const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
+                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
+                     for (index_t i=0; i < numComp; ++i) {
+                         o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
+                         o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
+                     } /* end of component loop i */
+                 } /* end of k0 loop */
+             } /* end of face 2 */
+             if (m_faceOffset[3] > -1) {
+ #pragma omp for nowait
+                 for (index_t k0=0; k0 < m_NE0; ++k0) {
+                     const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
+                     const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
+                     const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
+                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
+                     for (index_t i=0; i < numComp; ++i) {
+                         o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
+                         o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
+                     } /* end of component loop i */
+                 } /* end of k0 loop */
+             } /* end of face 3 */
+         } // end of parallel section
+     }
+ }
+ //protected
+ void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
+ {
+     const dim_t numComp = arg.getDataPointSize();
+     const double h0 = m_l0/m_gNE0;
+     const double h1 = m_l1/m_gNE1;
+     const index_t left = (m_offset0==0 ? 0 : 1);
+     const index_t bottom = (m_offset1==0 ? 0 : 1);
+     if (arg.getFunctionSpace().getTypeCode() == Elements) {
+         const double w = h0*h1/4.;
+ #pragma omp parallel
+         {
+             vector<double> int_local(numComp, 0);
+ #pragma omp for nowait
+             for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                     const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
+                     for (index_t i=0; i < numComp; ++i) {
+                         const double f0 = f[INDEX2(i,0,numComp)];
+                         const double f1 = f[INDEX2(i,1,numComp)];
+                         const double f2 = f[INDEX2(i,2,numComp)];
+                         const double f3 = f[INDEX2(i,3,numComp)];
+                         int_local[i]+=(f0+f1+f2+f3)*w;
+                     }  /* end of component loop i */
+                 } /* end of k0 loop */
+             } /* end of k1 loop */
+ #pragma omp critical
+             for (index_t i=0; i<numComp; i++)
+                 integrals[i]+=int_local[i];
+         } // end of parallel section
+     } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
+         const double w = h0*h1;
+ #pragma omp parallel
+         {
+             vector<double> int_local(numComp, 0);
+ #pragma omp for nowait
+             for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                     const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
+                     for (index_t i=0; i < numComp; ++i) {
+                         int_local[i]+=f[i]*w;
+                     }  /* end of component loop i */
+                 } /* end of k0 loop */
+             } /* end of k1 loop */
+ #pragma omp critical
+             for (index_t i=0; i<numComp; i++)
+                 integrals[i]+=int_local[i];
+         } // end of parallel section
+     } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
+         const double w0 = h0/2.;
+         const double w1 = h1/2.;
+ #pragma omp parallel
+         {
+             vector<double> int_local(numComp, 0);
+             if (m_faceOffset[0] > -1) {
+ #pragma omp for nowait
+                 for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                     const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
+                     for (index_t i=0; i < numComp; ++i) {
+                         const double f0 = f[INDEX2(i,0,numComp)];
+                         const double f1 = f[INDEX2(i,1,numComp)];
+                         int_local[i]+=(f0+f1)*w1;
+                     }  /* end of component loop i */
+                 } /* end of k1 loop */
+             }
+             if (m_faceOffset[1] > -1) {
+ #pragma omp for nowait
+                 for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                     const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
+                     for (index_t i=0; i < numComp; ++i) {
+                         const double f0 = f[INDEX2(i,0,numComp)];
+                         const double f1 = f[INDEX2(i,1,numComp)];
+                         int_local[i]+=(f0+f1)*w1;
+                     }  /* end of component loop i */
+                 } /* end of k1 loop */
+             }
+             if (m_faceOffset[2] > -1) {
+ #pragma omp for nowait
+                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                     const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
+                     for (index_t i=0; i < numComp; ++i) {
+                         const double f0 = f[INDEX2(i,0,numComp)];
+                         const double f1 = f[INDEX2(i,1,numComp)];
+                         int_local[i]+=(f0+f1)*w0;
+                     }  /* end of component loop i */
+                 } /* end of k0 loop */
+             }
+             if (m_faceOffset[3] > -1) {
+ #pragma omp for nowait
+                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                     const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
+                     for (index_t i=0; i < numComp; ++i) {
+                         const double f0 = f[INDEX2(i,0,numComp)];
+                         const double f1 = f[INDEX2(i,1,numComp)];
+                         int_local[i]+=(f0+f1)*w0;
+                     }  /* end of component loop i */
+                 } /* end of k0 loop */
+             }
+ #pragma omp critical
+             for (index_t i=0; i<numComp; i++)
+                 integrals[i]+=int_local[i];
+         } // end of parallel section
+     } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
+ #pragma omp parallel
+         {
+             vector<double> int_local(numComp, 0);
+             if (m_faceOffset[0] > -1) {
+ #pragma omp for nowait
+                 for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                     const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
+                     for (index_t i=0; i < numComp; ++i) {
+                         int_local[i]+=f[i]*h1;
+                     }  /* end of component loop i */
+                 } /* end of k1 loop */
+             }
+             if (m_faceOffset[1] > -1) {
+ #pragma omp for nowait
+                 for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
+                     const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
+                     for (index_t i=0; i < numComp; ++i) {
+                         int_local[i]+=f[i]*h1;
+                     }  /* end of component loop i */
+                 } /* end of k1 loop */
+             }
+             if (m_faceOffset[2] > -1) {
+ #pragma omp for nowait
+                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                     const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
+                     for (index_t i=0; i < numComp; ++i) {
+                         int_local[i]+=f[i]*h0;
+                     }  /* end of component loop i */
+                 } /* end of k0 loop */
+             }
+             if (m_faceOffset[3] > -1) {
+ #pragma omp for nowait
+                 for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
+                     const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
+                     for (index_t i=0; i < numComp; ++i) {
+                         int_local[i]+=f[i]*h0;
+                     }  /* end of component loop i */
+                 } /* end of k0 loop */
+             }
+ #pragma omp critical
+             for (index_t i=0; i<numComp; i++)
+                 integrals[i]+=int_local[i];
+         } // end of parallel section
+     }
+ }
+ //protected
  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
  {
      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
-Line 1110 
 void Rectangle::createPattern()
+Line 1111 
 void Rectangle::createPattern()
      // The rest is assigned in the loop further down
      m_dofMap.assign(getNumNodes(), 0);
  #pragma omp parallel for
-     for (index_t i=bottom; i<m_N1; i++) {
+     for (index_t i=bottom; i<bottom+nDOF1; i++) {
-         for (index_t j=left; j<m_N0; j++) {
+         for (index_t j=left; j<left+nDOF0; j++) {
              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
          }
      }
-Line 1272 
 void Rectangle::createPattern()
+Line 1273 
 void Rectangle::createPattern()
      Paso_Pattern_free(rowPattern);
  }
+ //private
+ void Rectangle::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
+          const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
+          bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
+ {
+     IndexVector rowIndex;
+     rowIndex.push_back(m_dofMap[firstNode]);
+     rowIndex.push_back(m_dofMap[firstNode+1]);
+     rowIndex.push_back(m_dofMap[firstNode+m_N0]);
+     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
+     if (addF) {
+         double *F_p=F.getSampleDataRW(0);
+         for (index_t i=0; i<rowIndex.size(); i++) {
+             if (rowIndex[i]<getNumDOF()) {
+                 for (index_t eq=0; eq<nEq; eq++) {
+                     F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
+                 }
+             }
+         }
+     }
+     if (addS) {
+         addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
+     }
+ }
  //protected
  void Rectangle::interpolateNodesOnElements(escript::Data& out,
                                          escript::Data& in, bool reduced) const
-Line 1279 
 void Rectangle::interpolateNodesOnElemen
+Line 1305 
 void Rectangle::interpolateNodesOnElemen
      const dim_t numComp = in.getDataPointSize();
      if (reduced) {
          out.requireWrite();
-         /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */
          const double c0 = .25;
  #pragma omp parallel for
          for (index_t k1=0; k1 < m_NE1; ++k1) {
              for (index_t k0=0; k0 < m_NE0; ++k0) {
-                 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));
-                 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));
-                 const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
+                 const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
-                 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));
                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
                  for (index_t i=0; i < numComp; ++i) {
                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
                  } /* end of component loop i */
              } /* end of k0 loop */
          } /* end of k1 loop */
-         /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */
      } else {
          out.requireWrite();
-         /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
          const double c0 = .16666666666666666667;
          const double c1 = .044658198738520451079;
          const double c2 = .62200846792814621559;
  #pragma omp parallel for
          for (index_t k1=0; k1 < m_NE1; ++k1) {
              for (index_t k0=0; k0 < m_NE0; ++k0) {
-                 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));
-                 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));
-                 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));
-                 const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
+                 const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
                  for (index_t i=0; i < numComp; ++i) {
                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);
-Line 1317 
 void Rectangle::interpolateNodesOnElemen
+Line 1340 
 void Rectangle::interpolateNodesOnElemen
                  } /* end of component loop i */
              } /* end of k0 loop */
          } /* end of k1 loop */
-         /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
      }
  }
-Line 1331 
 void Rectangle::interpolateNodesOnFaces(
+Line 1353 
 void Rectangle::interpolateNodesOnFaces(
          const double c0 = .5;
  #pragma omp parallel
          {
-             /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
                  for (index_t k1=0; k1 < m_NE1; ++k1) {
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
-                     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));
                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
-Line 1346 
 void Rectangle::interpolateNodesOnFaces(
+Line 1367 
 void Rectangle::interpolateNodesOnFaces(
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
                  for (index_t k1=0; k1 < m_NE1; ++k1) {
-                     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));
-                     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));
                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
-Line 1357 
 void Rectangle::interpolateNodesOnFaces(
+Line 1378 
 void Rectangle::interpolateNodesOnFaces(
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
                  for (index_t k0=0; k0 < m_NE0; ++k0) {
-                     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));
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
-Line 1368 
 void Rectangle::interpolateNodesOnFaces(
+Line 1389 
 void Rectangle::interpolateNodesOnFaces(
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
                  for (index_t k0=0; k0 < m_NE0; ++k0) {
-                     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));
-                     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));
                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
                      } /* end of component loop i */
                  } /* end of k0 loop */
              } /* end of face 3 */
-             /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */
          } // end of parallel section
      } else {
          out.requireWrite();
-Line 1384 
 void Rectangle::interpolateNodesOnFaces(
+Line 1404 
 void Rectangle::interpolateNodesOnFaces(
          const double c1 = 0.78867513459481288225;
  #pragma omp parallel
          {
-             /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */
              if (m_faceOffset[0] > -1) {
  #pragma omp for nowait
                  for (index_t k1=0; k1 < m_NE1; ++k1) {
-                     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));
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;
-Line 1400 
 void Rectangle::interpolateNodesOnFaces(
+Line 1419 
 void Rectangle::interpolateNodesOnFaces(
              if (m_faceOffset[1] > -1) {
  #pragma omp for nowait
                  for (index_t k1=0; k1 < m_NE1; ++k1) {
-                     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));
-                     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));
                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;
-Line 1412 
 void Rectangle::interpolateNodesOnFaces(
+Line 1431 
 void Rectangle::interpolateNodesOnFaces(
              if (m_faceOffset[2] > -1) {
  #pragma omp for nowait
                  for (index_t k0=0; k0 < m_NE0; ++k0) {
-                     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));
-                     const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
+                     const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;
-Line 1424 
 void Rectangle::interpolateNodesOnFaces(
+Line 1443 
 void Rectangle::interpolateNodesOnFaces(
              if (m_faceOffset[3] > -1) {
  #pragma omp for nowait
                  for (index_t k0=0; k0 < m_NE0; ++k0) {
-                     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));
-                     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));
                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
                      for (index_t i=0; i < numComp; ++i) {
                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;
-Line 1433 
 void Rectangle::interpolateNodesOnFaces(
+Line 1452 
 void Rectangle::interpolateNodesOnFaces(
                      } /* end of component loop i */
                  } /* end of k0 loop */
              } /* end of face 3 */
-             /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
          } // end of parallel section
      }
  }
-Line 1442 
 void Rectangle::interpolateNodesOnFaces(
+Line 1460 
 void Rectangle::interpolateNodesOnFaces(
  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
          const escript::Data& C, const escript::Data& D,
-         const escript::Data& X, const escript::Data& Y,
+         const escript::Data& X, const escript::Data& Y) const
-         const escript::Data& d, const escript::Data& y) const
  {
      const double h0 = m_l0/m_gNE0;
      const double h1 = m_l1/m_gNE1;
-     /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */
      const double w0 = -0.1555021169820365539*h1/h0;
      const double w1 = 0.041666666666666666667;
+     const double w2 = -0.15550211698203655390;
+     const double w3 = 0.041666666666666666667*h0/h1;
+     const double w4 = 0.15550211698203655390;
+     const double w5 = -0.041666666666666666667;
+     const double w6 = -0.01116454968463011277*h1/h0;
+     const double w7 = 0.011164549684630112770;
+     const double w8 = -0.011164549684630112770;
+     const double w9 = -0.041666666666666666667*h1/h0;
      const double w10 = -0.041666666666666666667*h0/h1;
      const double w11 = 0.1555021169820365539*h1/h0;
      const double w12 = 0.1555021169820365539*h0/h1;
-Line 1459 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1483 
 void Rectangle::assemblePDESingle(Paso_S
      const double w16 = -0.01116454968463011277*h0/h1;
      const double w17 = -0.1555021169820365539*h0/h1;
      const double w18 = -0.33333333333333333333*h1/h0;
-     const double w19 = 0.25000000000000000000;
+     const double w19 = 0.25;
-     const double w2 = -0.15550211698203655390;
+     const double w20 = -0.25;
-     const double w20 = -0.25000000000000000000;
      const double w21 = 0.16666666666666666667*h0/h1;
      const double w22 = -0.16666666666666666667*h1/h0;
      const double w23 = -0.16666666666666666667*h0/h1;
-Line 1471 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1494 
 void Rectangle::assemblePDESingle(Paso_S
      const double w27 = -0.33333333333333333333*h0/h1;
      const double w28 = -0.032861463941450536761*h1;
      const double w29 = -0.032861463941450536761*h0;
-     const double w3 = 0.041666666666666666667*h0/h1;
      const double w30 = -0.12264065304058601714*h1;
      const double w31 = -0.0023593469594139828636*h1;
      const double w32 = -0.008805202725216129906*h0;
-Line 1482 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1504 
 void Rectangle::assemblePDESingle(Paso_S
      const double w37 = 0.0023593469594139828636*h1;
      const double w38 = 0.12264065304058601714*h1;
      const double w39 = 0.032861463941450536761*h0;
-     const double w4 = 0.15550211698203655390;
      const double w40 = -0.12264065304058601714*h0;
      const double w41 = -0.0023593469594139828636*h0;
      const double w42 = 0.0023593469594139828636*h0;
-Line 1493 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1514 
 void Rectangle::assemblePDESingle(Paso_S
      const double w47 = 0.16666666666666666667*h1;
      const double w48 = 0.083333333333333333333*h0;
      const double w49 = -0.16666666666666666667*h0;
-     const double w5 = -0.041666666666666666667;
      const double w50 = 0.16666666666666666667*h0;
      const double w51 = -0.083333333333333333333*h1;
      const double w52 = 0.025917019497006092316*h0*h1;
-Line 1504 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1524 
 void Rectangle::assemblePDESingle(Paso_S
      const double w57 = 0.055555555555555555556*h0*h1;
      const double w58 = 0.027777777777777777778*h0*h1;
      const double w59 = 0.11111111111111111111*h0*h1;
-     const double w6 = -0.01116454968463011277*h1/h0;
      const double w60 = -0.19716878364870322056*h1;
      const double w61 = -0.19716878364870322056*h0;
      const double w62 = -0.052831216351296779436*h0;
-Line 1515 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1534 
 void Rectangle::assemblePDESingle(Paso_S
      const double w67 = 0.052831216351296779436*h0;
      const double w68 = -0.5*h1;
      const double w69 = -0.5*h0;
-     const double w7 = 0.011164549684630112770;
      const double w70 = 0.5*h1;
      const double w71 = 0.5*h0;
      const double w72 = 0.1555021169820365539*h0*h1;
      const double w73 = 0.041666666666666666667*h0*h1;
      const double w74 = 0.01116454968463011277*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 */
      rhs.requireWrite();
  #pragma omp parallel
-Line 1538 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1553 
 void Rectangle::assemblePDESingle(Paso_S
                      vector<double> EM_S(4*4, 0);
                      vector<double> EM_F(4, 0);
                      const index_t e = k0 + m_NE0*k1;
-                     /*** GENERATOR SNIP_PDE_SINGLE TOP */
                      ///////////////
                      // process A //
                      ///////////////
-Line 1546 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1560 
 void Rectangle::assemblePDESingle(Paso_S
                          add_EM_S=true;
                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
                          if (A.actsExpanded()) {
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             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)];
-                             const register double tmp4_0 = A_10_1 + A_10_2;
+                             const double tmp0_0 = A_01_0 + A_01_3;
-                             const register double tmp12_0 = A_11_0 + A_11_2;
+                             const double tmp1_0 = A_00_0 + A_00_1;
-                             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;
-                             const register double tmp10_0 = A_01_3 + A_10_3;
+                             const double tmp3_0 = A_00_2 + A_00_3;
-                             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;
-                             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;
-                             const register double tmp13_0 = A_01_2 + A_10_1;
+                             const double tmp6_0 = A_01_3 + A_10_0;
-                             const register double tmp3_0 = A_00_2 + A_00_3;
+                             const double tmp7_0 = A_01_0 + A_10_3;
-                             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;
-                             const register double tmp18_0 = A_01_1 + A_10_1;
+                             const double tmp9_0 = A_01_0 + A_10_0;
-                             const register double tmp1_0 = A_00_0 + A_00_1;
+                             const double tmp12_0 = A_11_0 + A_11_2;
-                             const register double tmp15_0 = A_01_1 + A_10_2;
+                             const double tmp10_0 = A_01_3 + A_10_3;
-                             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;
-                             const register double tmp16_0 = A_10_0 + A_10_3;
+                             const double tmp13_0 = A_01_2 + A_10_1;
-                             const register double tmp6_0 = A_01_3 + A_10_0;
+                             const double tmp11_0 = A_11_1 + A_11_3;
-                             const register double tmp17_0 = A_01_1 + A_01_2;
+                             const double tmp18_0 = A_01_1 + A_10_1;
-                             const register double tmp9_0 = A_01_0 + A_10_0;
+                             const double tmp15_0 = A_01_1 + A_10_2;
-                             const register double tmp7_0 = A_01_0 + A_10_3;
+                             const double tmp16_0 = A_10_0 + A_10_3;
-                             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;
-                             const register double tmp19_0 = A_01_2 + A_10_2;
+                             const double tmp19_0 = A_01_2 + A_10_2;
-                             const register double tmp14_1 = A_10_0*w8;
+                             const double tmp0_1 = A_10_3*w8;
-                             const register double tmp23_1 = tmp3_0*w14;
+                             const double tmp1_1 = tmp0_0*w1;
-                             const register double tmp35_1 = A_01_0*w8;
+                             const double tmp2_1 = A_01_1*w4;
-                             const register double tmp54_1 = tmp13_0*w8;
+                             const double tmp3_1 = tmp1_0*w0;
-                             const register double tmp20_1 = tmp9_0*w4;
+                             const double tmp4_1 = A_01_2*w7;
-                             const register double tmp25_1 = tmp12_0*w12;
+                             const double tmp5_1 = tmp2_0*w3;
-                             const register double tmp2_1 = A_01_1*w4;
+                             const double tmp6_1 = tmp3_0*w6;
-                             const register double tmp44_1 = tmp7_0*w7;
+                             const double tmp7_1 = A_10_0*w2;
-                             const register double tmp26_1 = tmp10_0*w4;
+                             const double tmp8_1 = tmp4_0*w5;
-                             const register double tmp52_1 = tmp18_0*w8;
+                             const double tmp9_1 = tmp2_0*w10;
-                             const register double tmp48_1 = A_10_1*w7;
+                             const double tmp14_1 = A_10_0*w8;
-                             const register double tmp46_1 = A_01_3*w8;
+                             const double tmp23_1 = tmp3_0*w14;
-                             const register double tmp50_1 = A_01_0*w2;
+                             const double tmp35_1 = A_01_0*w8;
-                             const register double tmp8_1 = tmp4_0*w5;
+                             const double tmp54_1 = tmp13_0*w8;
-                             const register double tmp56_1 = tmp19_0*w8;
+                             const double tmp20_1 = tmp9_0*w4;
-                             const register double tmp9_1 = tmp2_0*w10;
+                             const double tmp25_1 = tmp12_0*w12;
-                             const register double tmp19_1 = A_10_3*w2;
+                             const double tmp44_1 = tmp7_0*w7;
-                             const register double tmp47_1 = A_10_2*w4;
+                             const double tmp26_1 = tmp10_0*w4;
-                             const register double tmp16_1 = tmp3_0*w0;
+                             const double tmp52_1 = tmp18_0*w8;
-                             const register double tmp18_1 = tmp1_0*w6;
+                             const double tmp48_1 = A_10_1*w7;
-                             const register double tmp31_1 = tmp11_0*w12;
+                             const double tmp46_1 = A_01_3*w8;
-                             const register double tmp55_1 = tmp15_0*w2;
+                             const double tmp50_1 = A_01_0*w2;
-                             const register double tmp39_1 = A_10_2*w7;
+                             const double tmp56_1 = tmp19_0*w8;
-                             const register double tmp11_1 = tmp6_0*w7;
+                             const double tmp19_1 = A_10_3*w2;
-                             const register double tmp40_1 = tmp11_0*w17;
+                             const double tmp47_1 = A_10_2*w4;
-                             const register double tmp34_1 = tmp15_0*w8;
+                             const double tmp16_1 = tmp3_0*w0;
-                             const register double tmp33_1 = tmp14_0*w5;
+                             const double tmp18_1 = tmp1_0*w6;
-                             const register double tmp24_1 = tmp11_0*w13;
+                             const double tmp31_1 = tmp11_0*w12;
-                             const register double tmp3_1 = tmp1_0*w0;
+                             const double tmp55_1 = tmp15_0*w2;
-                             const register double tmp5_1 = tmp2_0*w3;
+                             const double tmp39_1 = A_10_2*w7;
-                             const register double tmp43_1 = tmp17_0*w5;
+                             const double tmp11_1 = tmp6_0*w7;
-                             const register double tmp15_1 = A_01_2*w4;
+                             const double tmp40_1 = tmp11_0*w17;
-                             const register double tmp53_1 = tmp19_0*w2;
+                             const double tmp34_1 = tmp15_0*w8;
-                             const register double tmp27_1 = tmp3_0*w11;
+                             const double tmp33_1 = tmp14_0*w5;
-                             const register double tmp32_1 = tmp13_0*w2;
+                             const double tmp24_1 = tmp11_0*w13;
-                             const register double tmp10_1 = tmp5_0*w9;
+                             const double tmp43_1 = tmp17_0*w5;
-                             const register double tmp37_1 = A_10_1*w4;
+                             const double tmp15_1 = A_01_2*w4;
-                             const register double tmp38_1 = tmp5_0*w15;
+                             const double tmp53_1 = tmp19_0*w2;
-                             const register double tmp17_1 = A_01_1*w7;
+                             const double tmp27_1 = tmp3_0*w11;
-                             const register double tmp12_1 = tmp7_0*w4;
+                             const double tmp32_1 = tmp13_0*w2;
-                             const register double tmp22_1 = tmp10_0*w7;
+                             const double tmp10_1 = tmp5_0*w9;
-                             const register double tmp57_1 = tmp18_0*w2;
+                             const double tmp37_1 = A_10_1*w4;
-                             const register double tmp28_1 = tmp9_0*w7;
+                             const double tmp38_1 = tmp5_0*w15;
-                             const register double tmp29_1 = tmp1_0*w14;
+                             const double tmp17_1 = A_01_1*w7;
-                             const register double tmp51_1 = tmp11_0*w16;
+                             const double tmp12_1 = tmp7_0*w4;
-                             const register double tmp42_1 = tmp12_0*w16;
+                             const double tmp22_1 = tmp10_0*w7;
-                             const register double tmp49_1 = tmp12_0*w17;
+                             const double tmp57_1 = tmp18_0*w2;
-                             const register double tmp21_1 = tmp1_0*w11;
+                             const double tmp28_1 = tmp9_0*w7;
-                             const register double tmp1_1 = tmp0_0*w1;
+                             const double tmp29_1 = tmp1_0*w14;
-                             const register double tmp45_1 = tmp6_0*w4;
+                             const double tmp51_1 = tmp11_0*w16;
-                             const register double tmp7_1 = A_10_0*w2;
+                             const double tmp42_1 = tmp12_0*w16;
-                             const register double tmp6_1 = tmp3_0*w6;
+                             const double tmp49_1 = tmp12_0*w17;
-                             const register double tmp13_1 = tmp8_0*w1;
+                             const double tmp21_1 = tmp1_0*w11;
-                             const register double tmp36_1 = tmp16_0*w1;
+                             const double tmp45_1 = tmp6_0*w4;
-                             const register double tmp41_1 = A_01_3*w2;
+                             const double tmp13_1 = tmp8_0*w1;
-                             const register double tmp30_1 = tmp12_0*w13;
+                             const double tmp36_1 = tmp16_0*w1;
-                             const register double tmp4_1 = A_01_2*w7;
+                             const double tmp41_1 = A_01_3*w2;
-                             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;
                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
-                             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;
+                             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(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
-                             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;
+                             EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
                              EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
+                             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,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
-                             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;
                              EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_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;
+                             EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
                              EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
-                             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;
-                         } 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;
-                             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;
-                             const register double A_01 = A_p[INDEX2(0,1,2)];
+                         } else { // constant data
-                             const register double A_10 = A_p[INDEX2(1,0,2)];
+                             const double A_00 = A_p[INDEX2(0,0,2)];
-                             const register double A_11 = A_p[INDEX2(1,1,2)];
+                             const double A_10 = A_p[INDEX2(1,0,2)];
-                             const register double tmp0_0 = A_01 + A_10;
+                             const double A_01 = A_p[INDEX2(0,1,2)];
-                             const register double tmp0_1 = A_00*w18;
+                             const double A_11 = A_p[INDEX2(1,1,2)];
-                             const register double tmp10_1 = A_01*w20;
+                             const double tmp0_0 = A_01 + A_10;
-                             const register double tmp12_1 = A_00*w26;
+                             const double tmp0_1 = A_00*w18;
-                             const register double tmp4_1 = A_00*w22;
+                             const double tmp1_1 = A_01*w19;
-                             const register double tmp8_1 = A_00*w24;
+                             const double tmp2_1 = A_10*w20;
-                             const register double tmp13_1 = A_10*w19;
+                             const double tmp3_1 = A_11*w21;
-                             const register double tmp9_1 = tmp0_0*w20;
+                             const double tmp4_1 = A_00*w22;
-                             const register double tmp3_1 = A_11*w21;
+                             const double tmp5_1 = tmp0_0*w19;
-                             const register double tmp11_1 = A_11*w27;
+                             const double tmp6_1 = A_11*w23;
-                             const register double tmp1_1 = A_01*w19;
+                             const double tmp7_1 = A_11*w25;
-                             const register double tmp6_1 = A_11*w23;
+                             const double tmp8_1 = A_00*w24;
-                             const register double tmp7_1 = A_11*w25;
+                             const double tmp9_1 = tmp0_0*w20;
-                             const register double tmp2_1 = A_10*w20;
+                             const double tmp10_1 = A_01*w20;
-                             const register double tmp5_1 = tmp0_0*w19;
+                             const double tmp11_1 = A_11*w27;
-                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
+                             const double tmp12_1 = A_00*w26;
-                             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;
                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
-                             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;
+                             EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
-                             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;
+                             EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
+                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
                              EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
-                             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;
                              EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
-                             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;
                              EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
-                             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;
+                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
+                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
                          }
                      }
                      ///////////////
-Line 1701 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1715 
 void Rectangle::assemblePDESingle(Paso_S
                          add_EM_S=true;
                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
                          if (B.actsExpanded()) {
-                             const register double B_0_0 = B_p[INDEX2(0,0,2)];
+                             const double B_0_0 = B_p[INDEX2(0,0,2)];
-                             const register double B_1_0 = B_p[INDEX2(1,0,2)];
+                             const double B_1_0 = B_p[INDEX2(1,0,2)];
-                             const register double B_0_1 = B_p[INDEX2(0,1,2)];
+                             const double B_0_1 = B_p[INDEX2(0,1,2)];
-                             const register double B_1_1 = B_p[INDEX2(1,1,2)];
+                             const double B_1_1 = B_p[INDEX2(1,1,2)];
-                             const register double B_0_2 = B_p[INDEX2(0,2,2)];
+                             const double B_0_2 = B_p[INDEX2(0,2,2)];
-                             const register double B_1_2 = B_p[INDEX2(1,2,2)];
+                             const double B_1_2 = B_p[INDEX2(1,2,2)];
-                             const register double B_0_3 = B_p[INDEX2(0,3,2)];
+                             const double B_0_3 = B_p[INDEX2(0,3,2)];
-                             const register double B_1_3 = B_p[INDEX2(1,3,2)];
+                             const double B_1_3 = B_p[INDEX2(1,3,2)];
-                             const register double tmp3_0 = B_0_0 + B_0_2;
+                             const double tmp0_0 = B_1_0 + B_1_1;
-                             const register double tmp1_0 = B_1_2 + B_1_3;
+                             const double tmp1_0 = B_1_2 + B_1_3;
-                             const register double tmp2_0 = B_0_1 + B_0_3;
+                             const double tmp2_0 = B_0_1 + B_0_3;
-                             const register double tmp0_0 = B_1_0 + B_1_1;
+                             const double tmp3_0 = B_0_0 + B_0_2;
-                             const register double tmp63_1 = B_1_1*w42;
+                             const double tmp63_1 = B_1_1*w42;
-                             const register double tmp79_1 = B_1_1*w40;
+                             const double tmp79_1 = B_1_1*w40;
-                             const register double tmp37_1 = tmp3_0*w35;
+                             const double tmp37_1 = tmp3_0*w35;
-                             const register double tmp8_1 = tmp0_0*w32;
+                             const double tmp8_1 = tmp0_0*w32;
-                             const register double tmp71_1 = B_0_1*w34;
+                             const double tmp71_1 = B_0_1*w34;
-                             const register double tmp19_1 = B_0_3*w31;
+                             const double tmp19_1 = B_0_3*w31;
-                             const register double tmp15_1 = B_0_3*w34;
+                             const double tmp15_1 = B_0_3*w34;
-                             const register double tmp9_1 = tmp3_0*w34;
+                             const double tmp9_1 = tmp3_0*w34;
-                             const register double tmp35_1 = B_1_0*w36;
+                             const double tmp35_1 = B_1_0*w36;
-                             const register double tmp66_1 = B_0_3*w28;
+                             const double tmp66_1 = B_0_3*w28;
-                             const register double tmp28_1 = B_1_0*w42;
+                             const double tmp28_1 = B_1_0*w42;
-                             const register double tmp22_1 = B_1_0*w40;
+                             const double tmp22_1 = B_1_0*w40;
-                             const register double tmp16_1 = B_1_2*w29;
+                             const double tmp16_1 = B_1_2*w29;
-                             const register double tmp6_1 = tmp2_0*w35;
+                             const double tmp6_1 = tmp2_0*w35;
-                             const register double tmp55_1 = B_1_3*w40;
+                             const double tmp55_1 = B_1_3*w40;
-                             const register double tmp50_1 = B_1_3*w42;
+                             const double tmp50_1 = B_1_3*w42;
-                             const register double tmp7_1 = tmp1_0*w29;
+                             const double tmp7_1 = tmp1_0*w29;
-                             const register double tmp1_1 = tmp1_0*w32;
+                             const double tmp1_1 = tmp1_0*w32;
-                             const register double tmp57_1 = B_0_3*w30;
+                             const double tmp57_1 = B_0_3*w30;
-                             const register double tmp18_1 = B_1_1*w32;
+                             const double tmp18_1 = B_1_1*w32;
-                             const register double tmp53_1 = B_1_0*w41;
+                             const double tmp53_1 = B_1_0*w41;
-                             const register double tmp61_1 = B_1_3*w36;
+                             const double tmp61_1 = B_1_3*w36;
-                             const register double tmp27_1 = B_0_3*w38;
+                             const double tmp27_1 = B_0_3*w38;
-                             const register double tmp64_1 = B_0_2*w30;
+                             const double tmp64_1 = B_0_2*w30;
-                             const register double tmp76_1 = B_0_1*w38;
+                             const double tmp76_1 = B_0_1*w38;
-                             const register double tmp39_1 = tmp2_0*w34;
+                             const double tmp39_1 = tmp2_0*w34;
-                             const register double tmp62_1 = B_0_1*w31;
+                             const double tmp62_1 = B_0_1*w31;
-                             const register double tmp56_1 = B_0_0*w31;
+                             const double tmp56_1 = B_0_0*w31;
-                             const register double tmp49_1 = B_1_1*w36;
+                             const double tmp49_1 = B_1_1*w36;
-                             const register double tmp2_1 = B_0_2*w31;
+                             const double tmp2_1 = B_0_2*w31;
-                             const register double tmp23_1 = B_0_2*w33;
+                             const double tmp23_1 = B_0_2*w33;
-                             const register double tmp38_1 = B_1_1*w43;
+                             const double tmp38_1 = B_1_1*w43;
-                             const register double tmp74_1 = B_1_2*w41;
+                             const double tmp74_1 = B_1_2*w41;
-                             const register double tmp43_1 = B_1_1*w41;
+                             const double tmp43_1 = B_1_1*w41;
-                             const register double tmp58_1 = B_0_2*w28;
+                             const double tmp58_1 = B_0_2*w28;
-                             const register double tmp67_1 = B_0_0*w33;
+                             const double tmp67_1 = B_0_0*w33;
-                             const register double tmp33_1 = tmp0_0*w39;
+                             const double tmp33_1 = tmp0_0*w39;
-                             const register double tmp4_1 = B_0_0*w28;
+                             const double tmp4_1 = B_0_0*w28;
-                             const register double tmp20_1 = B_0_0*w30;
+                             const double tmp20_1 = B_0_0*w30;
-                             const register double tmp13_1 = B_0_2*w38;
+                             const double tmp13_1 = B_0_2*w38;
-                             const register double tmp65_1 = B_1_2*w43;
+                             const double tmp65_1 = B_1_2*w43;
-                             const register double tmp0_1 = tmp0_0*w29;
+                             const double tmp0_1 = tmp0_0*w29;
-                             const register double tmp41_1 = tmp3_0*w33;
+                             const double tmp41_1 = tmp3_0*w33;
-                             const register double tmp73_1 = B_0_2*w37;
+                             const double tmp73_1 = B_0_2*w37;
-                             const register double tmp69_1 = B_0_0*w38;
+                             const double tmp69_1 = B_0_0*w38;
-                             const register double tmp48_1 = B_1_2*w39;
+                             const double tmp48_1 = B_1_2*w39;
-                             const register double tmp59_1 = B_0_1*w33;
+                             const double tmp59_1 = B_0_1*w33;
-                             const register double tmp17_1 = B_1_3*w41;
+                             const double tmp17_1 = B_1_3*w41;
-                             const register double tmp5_1 = B_0_3*w33;
+                             const double tmp5_1 = B_0_3*w33;
-                             const register double tmp3_1 = B_0_1*w30;
+                             const double tmp3_1 = B_0_1*w30;
-                             const register double tmp21_1 = B_0_1*w28;
+                             const double tmp21_1 = B_0_1*w28;
-                             const register double tmp42_1 = B_1_0*w29;
+                             const double tmp42_1 = B_1_0*w29;
-                             const register double tmp54_1 = B_1_2*w32;
+                             const double tmp54_1 = B_1_2*w32;
-                             const register double tmp60_1 = B_1_0*w39;
+                             const double tmp60_1 = B_1_0*w39;
-                             const register double tmp32_1 = tmp1_0*w36;
+                             const double tmp32_1 = tmp1_0*w36;
-                             const register double tmp10_1 = B_0_1*w37;
+                             const double tmp10_1 = B_0_1*w37;
-                             const register double tmp14_1 = B_0_0*w35;
+                             const double tmp14_1 = B_0_0*w35;
-                             const register double tmp29_1 = B_0_1*w35;
+                             const double tmp29_1 = B_0_1*w35;
-                             const register double tmp26_1 = B_1_2*w36;
+                             const double tmp26_1 = B_1_2*w36;
-                             const register double tmp30_1 = B_1_3*w43;
+                             const double tmp30_1 = B_1_3*w43;
-                             const register double tmp70_1 = B_0_2*w35;
+                             const double tmp70_1 = B_0_2*w35;
-                             const register double tmp34_1 = B_1_3*w39;
+                             const double tmp34_1 = B_1_3*w39;
-                             const register double tmp51_1 = B_1_0*w43;
+                             const double tmp51_1 = B_1_0*w43;
-                             const register double tmp31_1 = B_0_2*w34;
+                             const double tmp31_1 = B_0_2*w34;
-                             const register double tmp45_1 = tmp3_0*w28;
+                             const double tmp45_1 = tmp3_0*w28;
-                             const register double tmp11_1 = tmp1_0*w39;
+                             const double tmp11_1 = tmp1_0*w39;
-                             const register double tmp52_1 = B_1_1*w29;
+                             const double tmp52_1 = B_1_1*w29;
-                             const register double tmp44_1 = B_1_3*w32;
+                             const double tmp44_1 = B_1_3*w32;
-                             const register double tmp25_1 = B_1_1*w39;
+                             const double tmp25_1 = B_1_1*w39;
-                             const register double tmp47_1 = tmp2_0*w33;
+                             const double tmp47_1 = tmp2_0*w33;
-                             const register double tmp72_1 = B_1_3*w29;
+                             const double tmp72_1 = B_1_3*w29;
-                             const register double tmp40_1 = tmp2_0*w28;
+                             const double tmp40_1 = tmp2_0*w28;
-                             const register double tmp46_1 = B_1_2*w40;
+                             const double tmp46_1 = B_1_2*w40;
-                             const register double tmp36_1 = B_1_2*w42;
+                             const double tmp36_1 = B_1_2*w42;
-                             const register double tmp24_1 = B_0_0*w37;
+                             const double tmp24_1 = B_0_0*w37;
-                             const register double tmp77_1 = B_0_3*w35;
+                             const double tmp77_1 = B_0_3*w35;
-                             const register double tmp68_1 = B_0_3*w37;
+                             const double tmp68_1 = B_0_3*w37;
-                             const register double tmp78_1 = B_0_0*w34;
+                             const double tmp78_1 = B_0_0*w34;
-                             const register double tmp12_1 = tmp0_0*w36;
+                             const double tmp12_1 = tmp0_0*w36;
-                             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;
                              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(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;
+                             EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
-                             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;
+                             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(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
+                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
-                             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;
                              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(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;
                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
-                             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;
-                         } else { /* constant data */
+                             EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
-                             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;
-                             const register double B_1 = B_p[1];
+                         } else { // constant data
-                             const register double tmp6_1 = B_1*w50;
+                             const double B_0 = B_p[0];
-                             const register double tmp1_1 = B_1*w45;
+                             const double B_1 = B_p[1];
-                             const register double tmp5_1 = B_1*w49;
+                             const double tmp0_1 = B_0*w44;
-                             const register double tmp4_1 = B_1*w48;
+                             const double tmp1_1 = B_1*w45;
-                             const register double tmp0_1 = B_0*w44;
+                             const double tmp2_1 = B_0*w46;
-                             const register double tmp2_1 = B_0*w46;
+                             const double tmp3_1 = B_0*w47;
-                             const register double tmp7_1 = B_0*w51;
+                             const double tmp4_1 = B_1*w48;
-                             const register double tmp3_1 = B_0*w47;
+                             const double tmp5_1 = B_1*w49;
-                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
+                             const double tmp6_1 = B_1*w50;
-                             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;
                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
-                             EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
+                             EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
+                             EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
-                             EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
+                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
+                             EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
+                             EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
                              EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
-                             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;
                              EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
-                             EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
+                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
                              EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
-                             EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
+                             EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
+                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
+                             EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
                          }
                      }
                      ///////////////
-Line 1845 
 void Rectangle::assemblePDESingle(Paso_S
+Line 1859 
 void Rectangle::assemblePDESingle(Paso_S
                          add_EM_S=true;
                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
                          if (C.actsExpanded()) {
-                             const register double C_0_0 = C_p[INDEX2(0,0,2)];
+                             const double C_0_0 = C_p[INDEX2(0,0,2)];
-                             const register double C_1_0 = C_p[INDEX2(1,0,2)];
+                             const double C_1_0 = C_p[INDEX2(1,0,2)];
-                             const register double C_0_1 = C_p[INDEX2(0,1,2)];
+                             const double C_0_1 = C_p[INDEX2(0,1,2)];
-                             const register double C_1_1 = C_p[INDEX2(1,1,2)];
+                             const double C_1_1 = C_p[INDEX2(1,1,2)];
-                             const register double C_0_2 = C_p[INDEX2(0,2,2)];
+                             const double C_0_2 = C_p[INDEX2(0,2,2)];
-                             const register double C_1_2 = C_p[INDEX2(1,2,2)];
+                             const double C_1_2 = C_p[INDEX2(1,2,2)];
-                             const register double C_0_3 = C_p[INDEX2(0,3,2)];
+                             const double C_0_3 = C_p[INDEX2(0,3,2)];
-                             const register double C_1_3 = C_p[INDEX2(1,3,2)];
+                             const double C_1_3 = C_p[INDEX2(1,3,2)];
-                             const register double tmp2_0 = C_0_1 + C_0_3;
+                             const double tmp0_0 = C_1_0 + C_1_1;
-                             const register double tmp1_0 = C_1_2 + C_1_3;
+                             const double tmp1_0 = C_1_2 + C_1_3;
-                             const register double tmp3_0 = C_0_0 + C_0_2;
+                             const double tmp2_0 = C_0_1 + C_0_3;
-                             const register double tmp0_0 = C_1_0 + C_1_1;
+                             const double tmp3_0 = C_0_0 + C_0_2;
-                             const register double tmp64_1 = C_0_2*w30;
+                             const double tmp64_1 = C_0_2*w30;
-                             const register double tmp14_1 = C_0_2*w28;
+                             const double tmp14_1 = C_0_2*w28;
-                             const register double tmp19_1 = C_0_3*w31;
+                             const double tmp19_1 = C_0_3*w31;
-                             const register double tmp22_1 = C_1_0*w40;
+                             const double tmp22_1 = C_1_0*w40;
-                             const register double tmp37_1 = tmp3_0*w35;
+                             const double tmp37_1 = tmp3_0*w35;
-                             const register double tmp29_1 = C_0_1*w35;
+                             const double tmp29_1 = C_0_1*w35;
-                             const register double tmp73_1 = C_0_2*w37;
+                             const double tmp73_1 = C_0_2*w37;
-                             const register double tmp74_1 = C_1_2*w41;
+                             const double tmp74_1 = C_1_2*w41;
-                             const register double tmp52_1 = C_1_3*w39;
+                             const double tmp52_1 = C_1_3*w39;
-                             const register double tmp25_1 = C_1_1*w39;
+                             const double tmp25_1 = C_1_1*w39;
-                             const register double tmp62_1 = C_0_1*w31;
+                             const double tmp62_1 = C_0_1*w31;
-                             const register double tmp79_1 = C_1_1*w40;
+                             const double tmp79_1 = C_1_1*w40;
-                             const register double tmp43_1 = C_1_1*w36;
+                             const double tmp43_1 = C_1_1*w36;
-                             const register double tmp27_1 = C_0_3*w38;
+                             const double tmp27_1 = C_0_3*w38;
-                             const register double tmp28_1 = C_1_0*w42;
+                             const double tmp28_1 = C_1_0*w42;
-                             const register double tmp63_1 = C_1_1*w42;
+                             const double tmp63_1 = C_1_1*w42;
-                             const register double tmp59_1 = C_0_3*w34;
+                             const double tmp59_1 = C_0_3*w34;
-                             const register double tmp72_1 = C_1_3*w29;
+                             const double tmp72_1 = C_1_3*w29;
-                             const register double tmp40_1 = tmp2_0*w35;
+                             const double tmp40_1 = tmp2_0*w35;
-                             const register double tmp13_1 = C_0_3*w30;
+                             const double tmp13_1 = C_0_3*w30;
-                             const register double tmp51_1 = C_1_2*w40;
+                             const double tmp51_1 = C_1_2*w40;
-                             const register double tmp54_1 = C_1_2*w42;
+                             const double tmp54_1 = C_1_2*w42;
-                             const register double tmp12_1 = C_0_0*w31;
+                             const double tmp12_1 = C_0_0*w31;
-                             const register double tmp2_1 = tmp1_0*w32;
+                             const double tmp2_1 = tmp1_0*w32;
-                             const register double tmp68_1 = C_0_2*w31;
+                             const double tmp68_1 = C_0_2*w31;
-                             const register double tmp75_1 = C_1_0*w32;
+                             const double tmp75_1 = C_1_0*w32;
-                             const register double tmp49_1 = C_1_1*w41;
+                             const double tmp49_1 = C_1_1*w41;
-                             const register double tmp4_1 = C_0_2*w35;
+                             const double tmp4_1 = C_0_2*w35;
-                             const register double tmp66_1 = C_0_3*w28;
+                             const double tmp66_1 = C_0_3*w28;
-                             const register double tmp56_1 = C_0_1*w37;
+                             const double tmp56_1 = C_0_1*w37;
-                             const register double tmp5_1 = C_0_1*w34;
+                             const double tmp5_1 = C_0_1*w34;
-                             const register double tmp38_1 = tmp2_0*w34;
+                             const double tmp38_1 = tmp2_0*w34;
-                             const register double tmp76_1 = C_0_1*w38;
+                             const double tmp76_1 = C_0_1*w38;
-                             const register double tmp21_1 = C_0_1*w28;
+                             const double tmp21_1 = C_0_1*w28;
-                             const register double tmp69_1 = C_0_1*w30;
+                             const double tmp69_1 = C_0_1*w30;
-                             const register double tmp53_1 = C_1_0*w36;
+                             const double tmp53_1 = C_1_0*w36;
-                             const register double tmp42_1 = C_1_2*w39;
+                             const double tmp42_1 = C_1_2*w39;
-                             const register double tmp32_1 = tmp1_0*w29;
+                             const double tmp32_1 = tmp1_0*w29;
-                             const register double tmp45_1 = C_1_0*w43;
+                             const double tmp45_1 = C_1_0*w43;
-                             const register double tmp33_1 = tmp0_0*w32;
+                             const double tmp33_1 = tmp0_0*w32;
-                             const register double tmp35_1 = C_1_0*w41;
+                             const double tmp35_1 = C_1_0*w41;
-                             const register double tmp26_1 = C_1_2*w36;
+                             const double tmp26_1 = C_1_2*w36;
-                             const register double tmp67_1 = C_0_0*w33;
+                             const double tmp67_1 = C_0_0*w33;
-                             const register double tmp31_1 = C_0_2*w34;
+                             const double tmp31_1 = C_0_2*w34;
-                             const register double tmp20_1 = C_0_0*w30;
+                             const double tmp20_1 = C_0_0*w30;
-                             const register double tmp70_1 = C_0_0*w28;
+                             const double tmp70_1 = C_0_0*w28;
-                             const register double tmp8_1 = tmp0_0*w39;
+                             const double tmp8_1 = tmp0_0*w39;
-                             const register double tmp30_1 = C_1_3*w43;
+                             const double tmp30_1 = C_1_3*w43;
-                             const register double tmp0_1 = tmp0_0*w29;
+                             const double tmp0_1 = tmp0_0*w29;
-                             const register double tmp17_1 = C_1_3*w41;
+                             const double tmp17_1 = C_1_3*w41;
-                             const register double tmp58_1 = C_0_0*w35;
+                             const double tmp58_1 = C_0_0*w35;
-                             const register double tmp9_1 = tmp3_0*w33;
+                             const double tmp9_1 = tmp3_0*w33;
-                             const register double tmp61_1 = C_1_3*w36;
+                             const double tmp61_1 = C_1_3*w36;
-                             const register double tmp41_1 = tmp3_0*w34;
+                             const double tmp41_1 = tmp3_0*w34;
-                             const register double tmp50_1 = C_1_3*w32;
+                             const double tmp50_1 = C_1_3*w32;
-                             const register double tmp18_1 = C_1_1*w32;
+                             const double tmp18_1 = C_1_1*w32;
-                             const register double tmp6_1 = tmp1_0*w36;
+                             const double tmp6_1 = tmp1_0*w36;
-                             const register double tmp3_1 = C_0_0*w38;
+                             const double tmp3_1 = C_0_0*w38;
-                             const register double tmp34_1 = C_1_1*w29;
+                             const double tmp34_1 = C_1_1*w29;
-                             const register double tmp77_1 = C_0_3*w35;
+                             const double tmp77_1 = C_0_3*w35;
-                             const register double tmp65_1 = C_1_2*w43;
+                             const double tmp65_1 = C_1_2*w43;
-                             const register double tmp71_1 = C_0_3*w33;
+                             const double tmp71_1 = C_0_3*w33;
-                             const register double tmp55_1 = C_1_1*w43;
+                             const double tmp55_1 = C_1_1*w43;
-                             const register double tmp46_1 = tmp3_0*w28;
+                             const double tmp46_1 = tmp3_0*w28;
-                             const register double tmp24_1 = C_0_0*w37;
+                             const double tmp24_1 = C_0_0*w37;
-                             const register double tmp10_1 = tmp1_0*w39;
+                             const double tmp10_1 = tmp1_0*w39;
-                             const register double tmp48_1 = C_1_0*w29;
+                             const double tmp48_1 = C_1_0*w29;
-                             const register double tmp15_1 = C_0_1*w33;
+                             const double tmp15_1 = C_0_1*w33;
-                             const register double tmp36_1 = C_1_2*w32;
+                             const double tmp36_1 = C_1_2*w32;
-                             const register double tmp60_1 = C_1_0*w39;
+                             const double tmp60_1 = C_1_0*w39;
-                             const register double tmp47_1 = tmp2_0*w33;
+                             const double tmp47_1 = tmp2_0*w33;
-                             const register double tmp16_1 = C_1_2*w29;
+                             const double tmp16_1 = C_1_2*w29;
-                             const register double tmp1_1 = C_0_3*w37;
+                             const double tmp1_1 = C_0_3*w37;
-                             const register double tmp7_1 = tmp2_0*w28;
+                             const double tmp7_1 = tmp2_0*w28;
-                             const register double tmp39_1 = C_1_3*w40;
+                             const double tmp39_1 = C_1_3*w40;
-                             const register double tmp44_1 = C_1_3*w42;
+                             const double tmp44_1 = C_1_3*w42;
-                             const register double tmp57_1 = C_0_2*w38;
+                             const double tmp57_1 = C_0_2*w38;
-                             const register double tmp78_1 = C_0_0*w34;
+                             const double tmp78_1 = C_0_0*w34;
-                             const register double tmp11_1 = tmp0_0*w36;
+                             const double tmp11_1 = tmp0_0*w36;
-                             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;
                              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(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;
+                             EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
-                             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;
+                             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(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
+                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
-                             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;
                              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(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;
                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
-                             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;
-                         } else { /* constant data */
+                             EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
-                             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;
-                             const register double C_1 = C_p[1];
+                         } else { // constant data
-                             const register double tmp1_1 = C_1*w45;
+                             const double C_0 = C_p[0];
-                             const register double tmp3_1 = C_0*w51;
+                             const double C_1 = C_p[1];
-                             const register double tmp4_1 = C_0*w44;
+                             const double tmp0_1 = C_0*w47;
-                             const register double tmp7_1 = C_0*w46;
+                             const double tmp1_1 = C_1*w45;
-                             const register double tmp5_1 = C_1*w49;
+                             const double tmp2_1 = C_1*w48;
-                             const register double tmp2_1 = C_1*w48;
+                             const double tmp3_1 = C_0*w51;
-                             const register double tmp0_1 = C_0*w47;
+                             const double tmp4_1 = C_0*w44;
-                             const register double tmp6_1 = C_1*w50;
+                             const double tmp5_1 = C_1*w49;
-                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
+                             const double tmp6_1 = C_1*w50;
-                             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;
                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
-                             EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
+                             EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
+                             EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
-                             EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
+                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
+                             EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
                              EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
+                             EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
                              EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
-                             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;
                              EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
-                             EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
+                             EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
                              EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
-                             EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
+                             EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
+                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
+                             EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
                          }
                      }
                      ///////////////
-Line 1989 
 void Rectangle::assemblePDESingle(Paso_S
+Line 2003 
 void Rectangle::assemblePDESingle(Paso_S
                          add_EM_S=true;
                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
                          if (D.actsExpanded()) {
-                             const register double D_0 = D_p[0];
+                             const double D_0 = D_p[0];
-                             const register double D_1 = D_p[1];
+                             const double D_1 = D_p[1];
-                             const register double D_2 = D_p[2];
+                             const double D_2 = D_p[2];
-                             const register double D_3 = D_p[3];
+                             const double D_3 = D_p[3];
-                             const register double tmp4_0 = D_1 + D_3;
+                             const double tmp0_0 = D_0 + D_1;
-                             const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;
+                             const double tmp1_0 = D_2 + D_3;
-                             const register double tmp5_0 = D_0 + D_2;
+                             const double tmp2_0 = D_0 + D_1 + D_2 + D_3;
-                             const register double tmp0_0 = D_0 + D_1;
+                             const double tmp3_0 = D_1 + D_2;
-                             const register double tmp6_0 = D_0 + D_3;
+                             const double tmp4_0 = D_1 + D_3;
-                             const register double tmp1_0 = D_2 + D_3;
+                             const double tmp5_0 = D_0 + D_2;
-                             const register double tmp3_0 = D_1 + D_2;
+                             const double tmp6_0 = D_0 + D_3;
-                             const register double tmp16_1 = D_1*w56;
+                             const double tmp0_1 = tmp0_0*w52;
-                             const register double tmp14_1 = tmp6_0*w54;
+                             const double tmp1_1 = tmp1_0*w53;
-                             const register double tmp8_1 = D_3*w55;
+                             const double tmp2_1 = tmp2_0*w54;
-                             const register double tmp2_1 = tmp2_0*w54;
+                             const double tmp3_1 = tmp1_0*w52;
-                             const register double tmp12_1 = tmp5_0*w52;
+                             const double tmp4_1 = tmp0_0*w53;
-                             const register double tmp4_1 = tmp0_0*w53;
+                             const double tmp5_1 = tmp3_0*w54;
-                             const register double tmp3_1 = tmp1_0*w52;
+                             const double tmp6_1 = D_0*w55;
-                             const register double tmp13_1 = tmp4_0*w53;
+                             const double tmp7_1 = D_3*w56;
-                             const register double tmp10_1 = tmp4_0*w52;
+                             const double tmp8_1 = D_3*w55;
-                             const register double tmp1_1 = tmp1_0*w53;
+                             const double tmp9_1 = D_0*w56;
-                             const register double tmp7_1 = D_3*w56;
+                             const double tmp10_1 = tmp4_0*w52;
-                             const register double tmp0_1 = tmp0_0*w52;
+                             const double tmp11_1 = tmp5_0*w53;
-                             const register double tmp11_1 = tmp5_0*w53;
+                             const double tmp12_1 = tmp5_0*w52;
-                             const register double tmp9_1 = D_0*w56;
+                             const double tmp13_1 = tmp4_0*w53;
-                             const register double tmp5_1 = tmp3_0*w54;
+                             const double tmp14_1 = tmp6_0*w54;
-                             const register double tmp18_1 = D_2*w56;
+                             const double tmp15_1 = D_2*w55;
-                             const register double tmp17_1 = D_1*w55;
+                             const double tmp16_1 = D_1*w56;
-                             const register double tmp6_1 = D_0*w55;
+                             const double tmp17_1 = D_1*w55;
-                             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;
                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
-                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
+                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
+                             EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
                              EM_S[INDEX2(3,0,4)]+=tmp2_1;
-                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
+                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
+                             EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
                              EM_S[INDEX2(2,1,4)]+=tmp2_1;
+                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
                              EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
-                             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;
                              EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
-                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
+                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
                              EM_S[INDEX2(0,3,4)]+=tmp2_1;
-                             EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
+                             EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
-                         } else { /* constant data */
+                             EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
-                             const register double D_0 = D_p[0];
+                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
-                             const register double tmp2_1 = D_0*w59;
+                         } else { // constant data
-                             const register double tmp1_1 = D_0*w58;
+                             const double tmp0_1 = D_p[0]*w57;
-                             const register double tmp0_1 = D_0*w57;
+                             const double tmp1_1 = D_p[0]*w58;
-                             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;
                              EM_S[INDEX2(0,0,4)]+=tmp2_1;
-                             EM_S[INDEX2(3,3,4)]+=tmp2_1;
+                             EM_S[INDEX2(1,0,4)]+=tmp0_1;
+                             EM_S[INDEX2(2,0,4)]+=tmp0_1;
                              EM_S[INDEX2(3,0,4)]+=tmp1_1;
-                             EM_S[INDEX2(3,1,4)]+=tmp0_1;
+                             EM_S[INDEX2(0,1,4)]+=tmp0_1;
+                             EM_S[INDEX2(1,1,4)]+=tmp2_1;
                              EM_S[INDEX2(2,1,4)]+=tmp1_1;
+                             EM_S[INDEX2(3,1,4)]+=tmp0_1;
                              EM_S[INDEX2(0,2,4)]+=tmp0_1;
-                             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;
                              EM_S[INDEX2(2,2,4)]+=tmp2_1;
-                             EM_S[INDEX2(1,0,4)]+=tmp0_1;
+                             EM_S[INDEX2(3,2,4)]+=tmp0_1;
                              EM_S[INDEX2(0,3,4)]+=tmp1_1;
-                             EM_S[INDEX2(1,1,4)]+=tmp2_1;
+                             EM_S[INDEX2(1,3,4)]+=tmp0_1;
+                             EM_S[INDEX2(2,3,4)]+=tmp0_1;
+                             EM_S[INDEX2(3,3,4)]+=tmp2_1;
                          }
                      }
                      ///////////////
-Line 2065 
 void Rectangle::assemblePDESingle(Paso_S
+Line 2078 
 void Rectangle::assemblePDESingle(Paso_S
                          add_EM_F=true;
                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
                          if (X.actsExpanded()) {
-                             const register double X_0_0 = X_p[INDEX2(0,0,2)];
+                             const double X_0_0 = X_p[INDEX2(0,0,2)];
-                             const register double X_1_0 = X_p[INDEX2(1,0,2)];
+                             const double X_1_0 = X_p[INDEX2(1,0,2)];
-                             const register double X_0_1 = X_p[INDEX2(0,1,2)];
+                             const double X_0_1 = X_p[INDEX2(0,1,2)];
-                             const register double X_1_1 = X_p[INDEX2(1,1,2)];
+                             const double X_1_1 = X_p[INDEX2(1,1,2)];
-                             const register double X_0_2 = X_p[INDEX2(0,2,2)];
+                             const double X_0_2 = X_p[INDEX2(0,2,2)];
-                             const register double X_1_2 = X_p[INDEX2(1,2,2)];
+                             const double X_1_2 = X_p[INDEX2(1,2,2)];
-                             const register double X_0_3 = X_p[INDEX2(0,3,2)];
+                             const double X_0_3 = X_p[INDEX2(0,3,2)];
-                             const register double X_1_3 = X_p[INDEX2(1,3,2)];
+                             const double X_1_3 = X_p[INDEX2(1,3,2)];
-                             const register double tmp1_0 = X_1_1 + X_1_3;
+                             const double tmp0_0 = X_0_2 + X_0_3;
-                             const register double tmp3_0 = X_0_0 + X_0_1;
+                             const double tmp1_0 = X_1_1 + X_1_3;
-                             const register double tmp2_0 = X_1_0 + X_1_2;
+                             const double tmp2_0 = X_1_0 + X_1_2;
-                             const register double tmp0_0 = X_0_2 + X_0_3;
+                             const double tmp3_0 = X_0_0 + X_0_1;
-                             const register double tmp8_1 = tmp2_0*w66;
+                             const double tmp0_1 = tmp0_0*w63;
-                             const register double tmp5_1 = tmp3_0*w64;
+                             const double tmp1_1 = tmp1_0*w62;
-                             const register double tmp14_1 = tmp0_0*w64;
+                             const double tmp2_1 = tmp2_0*w61;
-                             const register double tmp3_1 = tmp3_0*w60;
+                             const double tmp3_1 = tmp3_0*w60;
-                             const register double tmp9_1 = tmp3_0*w63;
+                             const double tmp4_1 = tmp0_0*w65;
-                             const register double tmp13_1 = tmp3_0*w65;
+                             const double tmp5_1 = tmp3_0*w64;
-                             const register double tmp12_1 = tmp1_0*w66;
+                             const double tmp6_1 = tmp2_0*w62;
-                             const register double tmp10_1 = tmp0_0*w60;
+                             const double tmp7_1 = tmp1_0*w61;
-                             const register double tmp2_1 = tmp2_0*w61;
+                             const double tmp8_1 = tmp2_0*w66;
-                             const register double tmp6_1 = tmp2_0*w62;
+                             const double tmp9_1 = tmp3_0*w63;
-                             const register double tmp4_1 = tmp0_0*w65;
+                             const double tmp10_1 = tmp0_0*w60;
-                             const register double tmp11_1 = tmp1_0*w67;
+                             const double tmp11_1 = tmp1_0*w67;
-                             const register double tmp1_1 = tmp1_0*w62;
+                             const double tmp12_1 = tmp1_0*w66;
-                             const register double tmp7_1 = tmp1_0*w61;
+                             const double tmp13_1 = tmp3_0*w65;
-                             const register double tmp0_1 = tmp0_0*w63;
+                             const double tmp14_1 = tmp0_0*w64;
-                             const register double tmp15_1 = tmp2_0*w67;
+                             const double tmp15_1 = tmp2_0*w67;
                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
-                         } else { /* constant data */
+                         } else { // constant data
-                             const register double X_0 = X_p[0];
+                             const double tmp0_1 = X_p[1]*w69;
-                             const register double X_1 = X_p[1];
+                             const double tmp1_1 = X_p[0]*w68;
-                             const register double tmp3_1 = X_1*w71;
+                             const double tmp2_1 = X_p[0]*w70;
-                             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;
                              EM_F[0]+=tmp0_1 + tmp1_1;
                              EM_F[1]+=tmp0_1 + tmp2_1;
                              EM_F[2]+=tmp1_1 + tmp3_1;
-Line 2117 
 void Rectangle::assemblePDESingle(Paso_S
+Line 2128 
 void Rectangle::assemblePDESingle(Paso_S
                          add_EM_F=true;
                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
                          if (Y.actsExpanded()) {
-                             const register double Y_0 = Y_p[0];
+                             const double Y_0 = Y_p[0];
-                             const register double Y_1 = Y_p[1];
+                             const double Y_1 = Y_p[1];
-                             const register double Y_2 = Y_p[2];
+                             const double Y_2 = Y_p[2];
-                             const register double Y_3 = Y_p[3];
+                             const double Y_3 = Y_p[3];
-                             const register double tmp0_0 = Y_1 + Y_2;
+                             const double tmp0_0 = Y_1 + Y_2;
-                             const register double tmp1_0 = Y_0 + Y_3;
+                             const double tmp1_0 = Y_0 + Y_3;
-                             const register double tmp9_1 = Y_0*w74;
+                             const double tmp0_1 = Y_0*w72;
-                             const register double tmp4_1 = tmp1_0*w73;
+                             const double tmp1_1 = tmp0_0*w73;
-                             const register double tmp5_1 = Y_2*w74;
+                             const double tmp2_1 = Y_3*w74;
-                             const register double tmp7_1 = Y_1*w74;
+                             const double tmp3_1 = Y_1*w72;
-                             const register double tmp6_1 = Y_2*w72;
+                             const double tmp4_1 = tmp1_0*w73;
-                             const register double tmp2_1 = Y_3*w74;
+                             const double tmp5_1 = Y_2*w74;
-                             const register double tmp8_1 = Y_3*w72;
+                             const double tmp6_1 = Y_2*w72;
-                             const register double tmp3_1 = Y_1*w72;
+                             const double tmp7_1 = Y_1*w74;
-                             const register double tmp0_1 = Y_0*w72;
+                             const double tmp8_1 = Y_3*w72;
-                             const register double tmp1_1 = tmp0_0*w73;
+                             const double tmp9_1 = Y_0*w74;
                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
-                         } else { /* constant data */
+                         } else { // constant data
-                             const register double Y_0 = Y_p[0];
+                             const double tmp0_1 = Y_p[0]*w75;
-                             const register double tmp0_1 = Y_0*w75;
                              EM_F[0]+=tmp0_1;
                              EM_F[1]+=tmp0_1;
                              EM_F[2]+=tmp0_1;
                              EM_F[3]+=tmp0_1;
                          }
                      }
-                     /* GENERATOR SNIP_PDE_SINGLE BOTTOM */
                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
                      const index_t firstNode=m_N0*k1+k0;
-                     IndexVector rowIndex;
+                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
-                     rowIndex.push_back(m_dofMap[firstNode]);
+                 } // end k0 loop
-                     rowIndex.push_back(m_dofMap[firstNode+1]);
+             } // end k1 loop
-                     rowIndex.push_back(m_dofMap[firstNode+m_N0]);
+         } // end of colouring
-                     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
+     } // end of parallel region
-                     if (add_EM_F) {
+ }
-                         //cout << "-- AddtoRHS -- " << endl;
-                         double *F_p=rhs.getSampleDataRW(0);
+ //protected
-                         for (index_t i=0; i<rowIndex.size(); i++) {
+ void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
-                             if (rowIndex[i]<getNumDOF()) {
+         escript::Data& rhs, const escript::Data& A, const escript::Data& B,
-                                 F_p[rowIndex[i]]+=EM_F[i];
+         const escript::Data& C, const escript::Data& D,
-                                 //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;
+         const escript::Data& X, const escript::Data& Y) const
-                             }
+ {
-                         }
+     const double h0 = m_l0/m_gNE0;
-                         //cout << "---"<<endl;
+     const double h1 = m_l1/m_gNE1;
+     const double w0 = -.25*h1/h0;
+     const double w1 = .25;
+     const double w2 = -.25;
+     const double w3 = .25*h0/h1;
+     const double w4 = -.25*h0/h1;
+     const double w5 = .25*h1/h0;
+     const double w6 = -.125*h1;
+     const double w7 = -.125*h0;
+     const double w8 = .125*h1;
+     const double w9 = .125*h0;
+     const double w10 = .0625*h0*h1;
+     const double w11 = -.5*h1;
+     const double w12 = -.5*h0;
+     const double w13 = .5*h1;
+     const double w14 = .5*h0;
+     const double w15 = .25*h0*h1;
+     rhs.requireWrite();
+ #pragma omp parallel
+     {
+         for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
+ #pragma omp for
+             for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
+                 for (index_t k0=0; k0<m_NE0; ++k0)  {
+                     bool add_EM_S=false;
+                     bool add_EM_F=false;
+                     vector<double> EM_S(4*4, 0);
+                     vector<double> EM_F(4, 0);
+                     const index_t e = k0 + m_NE0*k1;
+                     ///////////////
+                     // process A //
+                     ///////////////
+                     if (!A.isEmpty()) {
+                         add_EM_S=true;
+                         const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
+                         const double A_00 = A_p[INDEX2(0,0,2)];
+                         const double A_10 = A_p[INDEX2(1,0,2)];
+                         const double A_01 = A_p[INDEX2(0,1,2)];
+                         const double A_11 = A_p[INDEX2(1,1,2)];
+                         const double tmp0_0 = A_01 + A_10;
+                         const double tmp0_1 = A_11*w3;
+                         const double tmp1_1 = A_00*w0;
+                         const double tmp2_1 = A_01*w1;
+                         const double tmp3_1 = A_10*w2;
+                         const double tmp4_1 = tmp0_0*w1;
+                         const double tmp5_1 = A_11*w4;
+                         const double tmp6_1 = A_00*w5;
+                         const double tmp7_1 = tmp0_0*w2;
+                         const double tmp8_1 = A_10*w1;
+                         const double tmp9_1 = A_01*w2;
+                         EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
+                         EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
+                         EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
+                         EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
+                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
+                         EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
+                         EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
+                         EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
+                         EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
+                         EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
+                         EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
+                         EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
+                         EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
+                         EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
+                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
+                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
                      }
-                     if (add_EM_S) {
+                     ///////////////
-                         //cout << "-- AddtoSystem -- " << endl;
+                     // process B //
-                         addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
+                     ///////////////
+                     if (!B.isEmpty()) {
+                         add_EM_S=true;
+                         const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
+                         const double tmp2_1 = B_p[0]*w8;
+                         const double tmp0_1 = B_p[0]*w6;
+                         const double tmp3_1 = B_p[1]*w9;
+                         const double tmp1_1 = B_p[1]*w7;
+                         EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
+                         EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
+                         EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
+                         EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;
+                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
+                         EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;
+                         EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;
+                         EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;
+                         EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;
+                         EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
+                         EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;
+                         EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
+                         EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;
+                         EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;
+                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;
+                         EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
+                     }
+                     ///////////////
+                     // process C //
+                     ///////////////
+                     if (!C.isEmpty()) {
+                         add_EM_S=true;
+                         const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
+                         const double tmp1_1 = C_p[1]*w7;
+                         const double tmp0_1 = C_p[0]*w8;
+                         const double tmp3_1 = C_p[0]*w6;
+                         const double tmp2_1 = C_p[1]*w9;
+                         EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
+                         EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
+                         EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
+                         EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
+                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
+                         EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;
+                         EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;
+                         EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;
+                         EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;
+                         EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
+                         EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
+                         EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
+                         EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;
+                         EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;
+                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
+                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;
+                     }
+                     ///////////////
+                     // process D //
+                     ///////////////
+                     if (!D.isEmpty()) {
+                         add_EM_S=true;
+                         const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
+                         const double tmp0_1 = D_p[0]*w10;
+                         EM_S[INDEX2(0,0,4)]+=tmp0_1;
+                         EM_S[INDEX2(1,0,4)]+=tmp0_1;
+                         EM_S[INDEX2(2,0,4)]+=tmp0_1;
+                         EM_S[INDEX2(3,0,4)]+=tmp0_1;
+                         EM_S[INDEX2(0,1,4)]+=tmp0_1;
+                         EM_S[INDEX2(1,1,4)]+=tmp0_1;
+                         EM_S[INDEX2(2,1,4)]+=tmp0_1;
+                         EM_S[INDEX2(3,1,4)]+=tmp0_1;
+                         EM_S[INDEX2(0,2,4)]+=tmp0_1;
+                         EM_S[INDEX2(1,2,4)]+=tmp0_1;
+                         EM_S[INDEX2(2,2,4)]+=tmp0_1;
+                         EM_S[INDEX2(3,2,4)]+=tmp0_1;
+                         EM_S[INDEX2(0,3,4)]+=tmp0_1;
+                         EM_S[INDEX2(1,3,4)]+=tmp0_1;
+                         EM_S[INDEX2(2,3,4)]+=tmp0_1;
+                         EM_S[INDEX2(3,3,4)]+=tmp0_1;
+                     }
+                     ///////////////
+                     // process X //
+                     ///////////////
+                     if (!X.isEmpty()) {
+                         add_EM_F=true;
+                         const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
+                         const double tmp0_1 = X_p[0]*w11;
+                         const double tmp2_1 = X_p[0]*w13;
+                         const double tmp1_1 = X_p[1]*w12;
+                         const double tmp3_1 = X_p[1]*w14;
+                         EM_F[0]+=tmp0_1 + tmp1_1;
+                         EM_F[1]+=tmp1_1 + tmp2_1;
+                         EM_F[2]+=tmp0_1 + tmp3_1;
+                         EM_F[3]+=tmp2_1 + tmp3_1;
+                     }
+                     ///////////////
+                     // process Y //
+                     ///////////////
+                     if (!Y.isEmpty()) {
+                         add_EM_F=true;
+                         const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
+                         const double tmp0_1 = Y_p[0]*w15;
+                         EM_F[0]+=tmp0_1;
+                         EM_F[1]+=tmp0_1;
+                         EM_F[2]+=tmp0_1;
+                         EM_F[3]+=tmp0_1;
                      }
+                     // add to matrix (if add_EM_S) and RHS (if add_EM_F)
+                     const index_t firstNode=m_N0*k1+k0;
+                     addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
                  } // end k0 loop
              } // end k1 loop
-         } // end of coloring
+         } // end of colouring
      } // end of parallel region
  }
-Line 2181 
 void Rectangle::assemblePDESingle(Paso_S
+Line 2360 
 void Rectangle::assemblePDESingle(Paso_S
  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
          const escript::Data& C, const escript::Data& D,
-         const escript::Data& X, const escript::Data& Y,
+         const escript::Data& X, const escript::Data& Y) const
-         const escript::Data& d, const escript::Data& y) const
  {
      const double h0 = m_l0/m_gNE0;
      const double h1 = m_l1/m_gNE1;
-Line 2194 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2372 
 void Rectangle::assemblePDESystem(Paso_S
          numComp=mat->logical_col_block_size;
      }
-     /* GENERATOR SNIP_PDE_SYSTEM_PRE TOP */
      const double w0 = -0.1555021169820365539*h1/h0;
      const double w1 = 0.041666666666666666667;
+     const double w2 = -0.15550211698203655390;
+     const double w3 = 0.041666666666666666667*h0/h1;
+     const double w4 = 0.15550211698203655390;
+     const double w5 = -0.041666666666666666667;
+     const double w6 = -0.01116454968463011277*h1/h0;
+     const double w7 = 0.011164549684630112770;
+     const double w8 = -0.011164549684630112770;
+     const double w9 = -0.041666666666666666667*h1/h0;
      const double w10 = -0.041666666666666666667*h0/h1;
      const double w11 = 0.1555021169820365539*h1/h0;
      const double w12 = 0.1555021169820365539*h0/h1;
-Line 2207 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2392 
 void Rectangle::assemblePDESystem(Paso_S
      const double w17 = -0.1555021169820365539*h0/h1;
      const double w18 = -0.33333333333333333333*h1/h0;
      const double w19 = 0.25000000000000000000;
-     const double w2 = -0.15550211698203655390;
      const double w20 = -0.25000000000000000000;
      const double w21 = 0.16666666666666666667*h0/h1;
      const double w22 = -0.16666666666666666667*h1/h0;
-Line 2218 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2402 
 void Rectangle::assemblePDESystem(Paso_S
      const double w27 = -0.33333333333333333333*h0/h1;
      const double w28 = -0.032861463941450536761*h1;
      const double w29 = -0.032861463941450536761*h0;
-     const double w3 = 0.041666666666666666667*h0/h1;
      const double w30 = -0.12264065304058601714*h1;
      const double w31 = -0.0023593469594139828636*h1;
      const double w32 = -0.008805202725216129906*h0;
-Line 2229 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2412 
 void Rectangle::assemblePDESystem(Paso_S
      const double w37 = 0.0023593469594139828636*h1;
      const double w38 = 0.12264065304058601714*h1;
      const double w39 = 0.032861463941450536761*h0;
-     const double w4 = 0.15550211698203655390;
      const double w40 = -0.12264065304058601714*h0;
      const double w41 = -0.0023593469594139828636*h0;
      const double w42 = 0.0023593469594139828636*h0;
-Line 2240 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2422 
 void Rectangle::assemblePDESystem(Paso_S
      const double w47 = 0.16666666666666666667*h1;
      const double w48 = 0.083333333333333333333*h0;
      const double w49 = -0.16666666666666666667*h0;
-     const double w5 = -0.041666666666666666667;
      const double w50 = 0.16666666666666666667*h0;
      const double w51 = -0.083333333333333333333*h1;
      const double w52 = 0.025917019497006092316*h0*h1;
-Line 2251 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2432 
 void Rectangle::assemblePDESystem(Paso_S
      const double w57 = 0.055555555555555555556*h0*h1;
      const double w58 = 0.027777777777777777778*h0*h1;
      const double w59 = 0.11111111111111111111*h0*h1;
-     const double w6 = -0.01116454968463011277*h1/h0;
      const double w60 = -0.19716878364870322056*h1;
      const double w61 = -0.19716878364870322056*h0;
      const double w62 = -0.052831216351296779436*h0;
-Line 2262 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2442 
 void Rectangle::assemblePDESystem(Paso_S
      const double w67 = 0.052831216351296779436*h0;
      const double w68 = -0.5*h1;
      const double w69 = -0.5*h0;
-     const double w7 = 0.011164549684630112770;
      const double w70 = 0.5*h1;
      const double w71 = 0.5*h0;
      const double w72 = 0.1555021169820365539*h0*h1;
      const double w73 = 0.041666666666666666667*h0*h1;
      const double w74 = 0.01116454968463011277*h0*h1;
      const double w75 = 0.25*h0*h1;
-     const double w8 = -0.011164549684630112770;
-     const double w9 = -0.041666666666666666667*h1/h0;
-     /* GENERATOR SNIP_PDE_SYSTEM_PRE BOTTOM */
      rhs.requireWrite();
  #pragma omp parallel
-Line 2285 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2461 
 void Rectangle::assemblePDESystem(Paso_S
                      vector<double> EM_S(4*4*numEq*numComp, 0);
                      vector<double> EM_F(4*numEq, 0);
                      const index_t e = k0 + m_NE0*k1;
-                     /* GENERATOR SNIP_PDE_SYSTEM TOP */
                      ///////////////
                      // process A //
                      ///////////////
-Line 2295 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2470 
 void Rectangle::assemblePDESystem(Paso_S
                          if (A.actsExpanded()) {
                              for (index_t k=0; k<numEq; k++) {
                                  for (index_t m=0; m<numComp; m++) {
-                                     const register double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];
+                                     const double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];
-                                     const register double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];
+                                     const double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];
-                                     const register double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];
+                                     const double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];
-                                     const register double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];
+                                     const double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];
-                                     const register double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];
+                                     const double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];
-                                     const register double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];
+                                     const double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];
-                                     const register double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];
+                                     const double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];
-                                     const register double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];
+                                     const double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];
-                                     const register double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];
+                                     const double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];
-                                     const register double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];
+                                     const double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];
-                                     const register double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];
+                                     const double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];
-                                     const register double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];
+                                     const double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];
-                                     const register double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];
+                                     const double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];
-                                     const register double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];
+                                     const double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];
-                                     const register double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];
+                                     const double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];
-                                     const register double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];
+                                     const double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];
-                                     const register double tmp4_0 = A_10_1 + A_10_2;
+                                     const double tmp0_0 = A_01_0 + A_01_3;
-                                     const register double tmp12_0 = A_11_0 + A_11_2;
+                                     const double tmp1_0 = A_00_0 + A_00_1;
-                                     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;
-                                     const register double tmp10_0 = A_01_3 + A_10_3;
+                                     const double tmp3_0 = A_00_2 + A_00_3;
-                                     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;
-                                     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;
-                                     const register double tmp13_0 = A_01_2 + A_10_1;
+                                     const double tmp6_0 = A_01_3 + A_10_0;
-                                     const register double tmp3_0 = A_00_2 + A_00_3;
+                                     const double tmp7_0 = A_01_0 + A_10_3;
-                                     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;
-                                     const register double tmp18_0 = A_01_1 + A_10_1;
+                                     const double tmp9_0 = A_01_0 + A_10_0;
-                                     const register double tmp1_0 = A_00_0 + A_00_1;
+                                     const double tmp10_0 = A_01_3 + A_10_3;
-                                     const register double tmp15_0 = A_01_1 + A_10_2;
+                                     const double tmp11_0 = A_11_1 + A_11_3;
-                                     const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
+                                     const double tmp12_0 = A_11_0 + A_11_2;
-                                     const register double tmp16_0 = A_10_0 + A_10_3;
+                                     const double tmp13_0 = A_01_2 + A_10_1;
-                                     const register double tmp6_0 = A_01_3 + A_10_0;
+                                     const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
-                                     const register double tmp17_0 = A_01_1 + A_01_2;
+                                     const double tmp15_0 = A_01_1 + A_10_2;
-                                     const register double tmp9_0 = A_01_0 + A_10_0;
+                                     const double tmp16_0 = A_10_0 + A_10_3;
-                                     const register double tmp7_0 = A_01_0 + A_10_3;
+                                     const double tmp17_0 = A_01_1 + A_01_2;
-                                     const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
+                                     const double tmp18_0 = A_01_1 + A_10_1;
-                                     const register double tmp19_0 = A_01_2 + A_10_2;
+                                     const double tmp19_0 = A_01_2 + A_10_2;
-                                     const register double tmp14_1 = A_10_0*w8;
+                                     const double tmp0_1 = A_10_3*w8;
-                                     const register double tmp23_1 = tmp3_0*w14;
+                                     const double tmp1_1 = tmp0_0*w1;
-                                     const register double tmp35_1 = A_01_0*w8;
+                                     const double tmp2_1 = A_01_1*w4;
-                                     const register double tmp54_1 = tmp13_0*w8;
+                                     const double tmp3_1 = tmp1_0*w0;
-                                     const register double tmp20_1 = tmp9_0*w4;
+                                     const double tmp4_1 = A_01_2*w7;
-                                     const register double tmp25_1 = tmp12_0*w12;
+                                     const double tmp5_1 = tmp2_0*w3;
-                                     const register double tmp2_1 = A_01_1*w4;
+                                     const double tmp6_1 = tmp3_0*w6;
-                                     const register double tmp44_1 = tmp7_0*w7;
+                                     const double tmp7_1 = A_10_0*w2;
-                                     const register double tmp26_1 = tmp10_0*w4;
+                                     const double tmp8_1 = tmp4_0*w5;
-                                     const register double tmp52_1 = tmp18_0*w8;
+                                     const double tmp9_1 = tmp2_0*w10;
-                                     const register double tmp48_1 = A_10_1*w7;
+                                     const double tmp10_1 = tmp5_0*w9;
-                                     const register double tmp46_1 = A_01_3*w8;
+                                     const double tmp11_1 = tmp6_0*w7;
-                                     const register double tmp50_1 = A_01_0*w2;
+                                     const double tmp12_1 = tmp7_0*w4;
-                                     const register double tmp8_1 = tmp4_0*w5;
+                                     const double tmp13_1 = tmp8_0*w1;
-                                     const register double tmp56_1 = tmp19_0*w8;
+                                     const double tmp14_1 = A_10_0*w8;
-                                     const register double tmp9_1 = tmp2_0*w10;
+                                     const double tmp15_1 = A_01_2*w4;
-                                     const register double tmp19_1 = A_10_3*w2;
+                                     const double tmp16_1 = tmp3_0*w0;
-                                     const register double tmp47_1 = A_10_2*w4;
+                                     const double tmp17_1 = A_01_1*w7;
-                                     const register double tmp16_1 = tmp3_0*w0;
+                                     const double tmp18_1 = tmp1_0*w6;
-                                     const register double tmp18_1 = tmp1_0*w6;
+                                     const double tmp19_1 = A_10_3*w2;
-                                     const register double tmp31_1 = tmp11_0*w12;
+                                     const double tmp20_1 = tmp9_0*w4;
-                                     const register double tmp55_1 = tmp15_0*w2;
+                                     const double tmp21_1 = tmp1_0*w11;
-                                     const register double tmp39_1 = A_10_2*w7;
+                                     const double tmp22_1 = tmp10_0*w7;
-                                     const register double tmp11_1 = tmp6_0*w7;
+                                     const double tmp23_1 = tmp3_0*w14;
-                                     const register double tmp40_1 = tmp11_0*w17;
+                                     const double tmp24_1 = tmp11_0*w13;
-                                     const register double tmp34_1 = tmp15_0*w8;
+                                     const double tmp25_1 = tmp12_0*w12;
-                                     const register double tmp33_1 = tmp14_0*w5;
+                                     const double tmp26_1 = tmp10_0*w4;
-                                     const register double tmp24_1 = tmp11_0*w13;
+                                     const double tmp27_1 = tmp3_0*w11;
-                                     const register double tmp3_1 = tmp1_0*w0;
+                                     const double tmp28_1 = tmp9_0*w7;
-                                     const register double tmp5_1 = tmp2_0*w3;
+                                     const double tmp29_1 = tmp1_0*w14;
-                                     const register double tmp43_1 = tmp17_0*w5;
+                                     const double tmp30_1 = tmp12_0*w13;
-                                     const register double tmp15_1 = A_01_2*w4;
+                                     const double tmp31_1 = tmp11_0*w12;
-                                     const register double tmp53_1 = tmp19_0*w2;
+                                     const double tmp32_1 = tmp13_0*w2;
-                                     const register double tmp27_1 = tmp3_0*w11;
+                                     const double tmp33_1 = tmp14_0*w5;
-                                     const register double tmp32_1 = tmp13_0*w2;
+                                     const double tmp34_1 = tmp15_0*w8;
-                                     const register double tmp10_1 = tmp5_0*w9;
+                                     const double tmp35_1 = A_01_0*w8;
-                                     const register double tmp37_1 = A_10_1*w4;
+                                     const double tmp36_1 = tmp16_0*w1;
-                                     const register double tmp38_1 = tmp5_0*w15;
+                                     const double tmp37_1 = A_10_1*w4;
-                                     const register double tmp17_1 = A_01_1*w7;
+                                     const double tmp38_1 = tmp5_0*w15;
-                                     const register double tmp12_1 = tmp7_0*w4;
+                                     const double tmp39_1 = A_10_2*w7;
-                                     const register double tmp22_1 = tmp10_0*w7;
+                                     const double tmp40_1 = tmp11_0*w17;
-                                     const register double tmp57_1 = tmp18_0*w2;
+                                     const double tmp41_1 = A_01_3*w2;
-                                     const register double tmp28_1 = tmp9_0*w7;
+                                     const double tmp42_1 = tmp12_0*w16;
-                                     const register double tmp29_1 = tmp1_0*w14;
+                                     const double tmp43_1 = tmp17_0*w5;
-                                     const register double tmp51_1 = tmp11_0*w16;
+                                     const double tmp44_1 = tmp7_0*w7;
-                                     const register double tmp42_1 = tmp12_0*w16;
+                                     const double tmp45_1 = tmp6_0*w4;
-                                     const register double tmp49_1 = tmp12_0*w17;
+                                     const double tmp46_1 = A_01_3*w8;
-                                     const register double tmp21_1 = tmp1_0*w11;
+                                     const double tmp47_1 = A_10_2*w4;
-                                     const register double tmp1_1 = tmp0_0*w1;
+                                     const double tmp48_1 = A_10_1*w7;
-                                     const register double tmp45_1 = tmp6_0*w4;
+                                     const double tmp49_1 = tmp12_0*w17;
-                                     const register double tmp7_1 = A_10_0*w2;
+                                     const double tmp50_1 = A_01_0*w2;
-                                     const register double tmp6_1 = tmp3_0*w6;
+                                     const double tmp51_1 = tmp11_0*w16;
-                                     const register double tmp13_1 = tmp8_0*w1;
+                                     const double tmp52_1 = tmp18_0*w8;
-                                     const register double tmp36_1 = tmp16_0*w1;
+                                     const double tmp53_1 = tmp19_0*w2;
-                                     const register double tmp41_1 = A_01_3*w2;
+                                     const double tmp54_1 = tmp13_0*w8;
-                                     const register double tmp30_1 = tmp12_0*w13;
+                                     const double tmp55_1 = tmp15_0*w2;
-                                     const register double tmp4_1 = A_01_2*w7;
+                                     const double tmp56_1 = tmp19_0*w8;
-                                     const register double tmp0_1 = A_10_3*w8;
+                                     const double tmp57_1 = tmp18_0*w2;
-                                     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;
-                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
-                                     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;
                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
-                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
+                                     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;
+                                     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;
                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
-                                     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;
+                                     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;
+                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
+                                     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;
                                      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;
-                                     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;
+                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
-                                     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;
-                                     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;
                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
-                                     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;
+                                     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;
                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
-                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
+                                     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;
+                                     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;
+                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
                                  }
                              }
-                         } else { /* constant data */
+                         } else { // constant data
                              for (index_t k=0; k<numEq; k++) {
                                  for (index_t m=0; m<numComp; m++) {
-                                     const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
+                                     const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
-                                     const register double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
+                                     const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
-                                     const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
+                                     const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
-                                     const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
+                                     const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
-                                     const register double tmp0_0 = A_01 + A_10;
+                                     const double tmp0_0 = A_01 + A_10;
-                                     const register double tmp0_1 = A_00*w18;
+                                     const double tmp0_1 = A_00*w18;
-                                     const register double tmp10_1 = A_01*w20;
+                                     const double tmp1_1 = A_01*w19;
-                                     const register double tmp12_1 = A_00*w26;
+                                     const double tmp2_1 = A_10*w20;
-                                     const register double tmp4_1 = A_00*w22;
+                                     const double tmp3_1 = A_11*w21;
-                                     const register double tmp8_1 = A_00*w24;
+                                     const double tmp4_1 = A_00*w22;
-                                     const register double tmp13_1 = A_10*w19;
+                                     const double tmp5_1 = tmp0_0*w19;
-                                     const register double tmp9_1 = tmp0_0*w20;
+                                     const double tmp6_1 = A_11*w23;
-                                     const register double tmp3_1 = A_11*w21;
+                                     const double tmp7_1 = A_11*w25;
-                                     const register double tmp11_1 = A_11*w27;
+                                     const double tmp8_1 = A_00*w24;
-                                     const register double tmp1_1 = A_01*w19;
+                                     const double tmp9_1 = tmp0_0*w20;
-                                     const register double tmp6_1 = A_11*w23;
+                                     const double tmp10_1 = A_01*w20;
-                                     const register double tmp7_1 = A_11*w25;
+                                     const double tmp11_1 = A_11*w27;
-                                     const register double tmp2_1 = A_10*w20;
+                                     const double tmp12_1 = A_00*w26;
-                                     const register double tmp5_1 = tmp0_0*w19;
+                                     const double tmp13_1 = A_10*w19;
-                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
-                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
-                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
-                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
+                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
+                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
-                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
+                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
+                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
+                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
-                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
+                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
-                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
-                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
-                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
+                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
-                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
+                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
+                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
+                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
                                  }
                              }
                          }
-Line 2458 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2633 
 void Rectangle::assemblePDESystem(Paso_S
                          if (B.actsExpanded()) {
                              for (index_t k=0; k<numEq; k++) {
                                  for (index_t m=0; m<numComp; m++) {
-                                     const register double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
+                                     const double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
-                                     const register double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
+                                     const double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
-                                     const register double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
+                                     const double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
-                                     const register double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
+                                     const double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
-                                     const register double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];
+                                     const double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];
-                                     const register double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
+                                     const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
-                                     const register double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
+                                     const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
-                                     const register double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
+                                     const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
-                                     const register double tmp3_0 = B_0_0 + B_0_2;
+                                     const double tmp0_0 = B_1_0 + B_1_1;
-                                     const register double tmp1_0 = B_1_2 + B_1_3;
+                                     const double tmp1_0 = B_1_2 + B_1_3;
-                                     const register double tmp2_0 = B_0_1 + B_0_3;
+                                     const double tmp2_0 = B_0_1 + B_0_3;
-                                     const register double tmp0_0 = B_1_0 + B_1_1;
+                                     const double tmp3_0 = B_0_0 + B_0_2;
-                                     const register double tmp63_1 = B_1_1*w42;
+                                     const double tmp63_1 = B_1_1*w42;
-                                     const register double tmp79_1 = B_1_1*w40;
+                                     const double tmp79_1 = B_1_1*w40;
-                                     const register double tmp37_1 = tmp3_0*w35;
+                                     const double tmp37_1 = tmp3_0*w35;
-                                     const register double tmp8_1 = tmp0_0*w32;
+                                     const double tmp8_1 = tmp0_0*w32;
-                                     const register double tmp71_1 = B_0_1*w34;
+                                     const double tmp71_1 = B_0_1*w34;
-                                     const register double tmp19_1 = B_0_3*w31;
+                                     const double tmp19_1 = B_0_3*w31;
-                                     const register double tmp15_1 = B_0_3*w34;
+                                     const double tmp15_1 = B_0_3*w34;
-                                     const register double tmp9_1 = tmp3_0*w34;
+                                     const double tmp9_1 = tmp3_0*w34;
-                                     const register double tmp35_1 = B_1_0*w36;
+                                     const double tmp35_1 = B_1_0*w36;
-                                     const register double tmp66_1 = B_0_3*w28;
+                                     const double tmp66_1 = B_0_3*w28;
-                                     const register double tmp28_1 = B_1_0*w42;
+                                     const double tmp28_1 = B_1_0*w42;
-                                     const register double tmp22_1 = B_1_0*w40;
+                                     const double tmp22_1 = B_1_0*w40;
-                                     const register double tmp16_1 = B_1_2*w29;
+                                     const double tmp16_1 = B_1_2*w29;
-                                     const register double tmp6_1 = tmp2_0*w35;
+                                     const double tmp6_1 = tmp2_0*w35;
-                                     const register double tmp55_1 = B_1_3*w40;
+                                     const double tmp55_1 = B_1_3*w40;
-                                     const register double tmp50_1 = B_1_3*w42;
+                                     const double tmp50_1 = B_1_3*w42;
-                                     const register double tmp7_1 = tmp1_0*w29;
+                                     const double tmp7_1 = tmp1_0*w29;
-                                     const register double tmp1_1 = tmp1_0*w32;
+                                     const double tmp1_1 = tmp1_0*w32;
-                                     const register double tmp57_1 = B_0_3*w30;
+                                     const double tmp57_1 = B_0_3*w30;
-                                     const register double tmp18_1 = B_1_1*w32;
+                                     const double tmp18_1 = B_1_1*w32;
-                                     const register double tmp53_1 = B_1_0*w41;
+                                     const double tmp53_1 = B_1_0*w41;
-                                     const register double tmp61_1 = B_1_3*w36;
+                                     const double tmp61_1 = B_1_3*w36;
-                                     const register double tmp27_1 = B_0_3*w38;
+                                     const double tmp27_1 = B_0_3*w38;
-                                     const register double tmp64_1 = B_0_2*w30;
+                                     const double tmp64_1 = B_0_2*w30;
-                                     const register double tmp76_1 = B_0_1*w38;
+                                     const double tmp76_1 = B_0_1*w38;
-                                     const register double tmp39_1 = tmp2_0*w34;
+                                     const double tmp39_1 = tmp2_0*w34;
-                                     const register double tmp62_1 = B_0_1*w31;
+                                     const double tmp62_1 = B_0_1*w31;
-                                     const register double tmp56_1 = B_0_0*w31;
+                                     const double tmp56_1 = B_0_0*w31;
-                                     const register double tmp49_1 = B_1_1*w36;
+                                     const double tmp49_1 = B_1_1*w36;
-                                     const register double tmp2_1 = B_0_2*w31;
+                                     const double tmp2_1 = B_0_2*w31;
-                                     const register double tmp23_1 = B_0_2*w33;
+                                     const double tmp23_1 = B_0_2*w33;
-                                     const register double tmp38_1 = B_1_1*w43;
+                                     const double tmp38_1 = B_1_1*w43;
-                                     const register double tmp74_1 = B_1_2*w41;
+                                     const double tmp74_1 = B_1_2*w41;
-                                     const register double tmp43_1 = B_1_1*w41;
+                                     const double tmp43_1 = B_1_1*w41;
-                                     const register double tmp58_1 = B_0_2*w28;
+                                     const double tmp58_1 = B_0_2*w28;
-                                     const register double tmp67_1 = B_0_0*w33;
+                                     const double tmp67_1 = B_0_0*w33;
-                                     const register double tmp33_1 = tmp0_0*w39;
+                                     const double tmp33_1 = tmp0_0*w39;
-                                     const register double tmp4_1 = B_0_0*w28;
+                                     const double tmp4_1 = B_0_0*w28;
-                                     const register double tmp20_1 = B_0_0*w30;
+                                     const double tmp20_1 = B_0_0*w30;
-                                     const register double tmp13_1 = B_0_2*w38;
+                                     const double tmp13_1 = B_0_2*w38;
-                                     const register double tmp65_1 = B_1_2*w43;
+                                     const double tmp65_1 = B_1_2*w43;
-                                     const register double tmp0_1 = tmp0_0*w29;
+                                     const double tmp0_1 = tmp0_0*w29;
-                                     const register double tmp41_1 = tmp3_0*w33;
+                                     const double tmp41_1 = tmp3_0*w33;
-                                     const register double tmp73_1 = B_0_2*w37;
+                                     const double tmp73_1 = B_0_2*w37;
-                                     const register double tmp69_1 = B_0_0*w38;
+                                     const double tmp69_1 = B_0_0*w38;
-                                     const register double tmp48_1 = B_1_2*w39;
+                                     const double tmp48_1 = B_1_2*w39;
-                                     const register double tmp59_1 = B_0_1*w33;
+                                     const double tmp59_1 = B_0_1*w33;
-                                     const register double tmp17_1 = B_1_3*w41;
+                                     const double tmp17_1 = B_1_3*w41;
-                                     const register double tmp5_1 = B_0_3*w33;
+                                     const double tmp5_1 = B_0_3*w33;
-                                     const register double tmp3_1 = B_0_1*w30;
+                                     const double tmp3_1 = B_0_1*w30;
-                                     const register double tmp21_1 = B_0_1*w28;
+                                     const double tmp21_1 = B_0_1*w28;
-                                     const register double tmp42_1 = B_1_0*w29;
+                                     const double tmp42_1 = B_1_0*w29;
-                                     const register double tmp54_1 = B_1_2*w32;
+                                     const double tmp54_1 = B_1_2*w32;
-                                     const register double tmp60_1 = B_1_0*w39;
+                                     const double tmp60_1 = B_1_0*w39;
-                                     const register double tmp32_1 = tmp1_0*w36;
+                                     const double tmp32_1 = tmp1_0*w36;
-                                     const register double tmp10_1 = B_0_1*w37;
+                                     const double tmp10_1 = B_0_1*w37;
-                                     const register double tmp14_1 = B_0_0*w35;
+                                     const double tmp14_1 = B_0_0*w35;
-                                     const register double tmp29_1 = B_0_1*w35;
+                                     const double tmp29_1 = B_0_1*w35;
-                                     const register double tmp26_1 = B_1_2*w36;
+                                     const double tmp26_1 = B_1_2*w36;
-                                     const register double tmp30_1 = B_1_3*w43;
+                                     const double tmp30_1 = B_1_3*w43;
-                                     const register double tmp70_1 = B_0_2*w35;
+                                     const double tmp70_1 = B_0_2*w35;
-                                     const register double tmp34_1 = B_1_3*w39;
+                                     const double tmp34_1 = B_1_3*w39;
-                                     const register double tmp51_1 = B_1_0*w43;
+                                     const double tmp51_1 = B_1_0*w43;
-                                     const register double tmp31_1 = B_0_2*w34;
+                                     const double tmp31_1 = B_0_2*w34;
-                                     const register double tmp45_1 = tmp3_0*w28;
+                                     const double tmp45_1 = tmp3_0*w28;
-                                     const register double tmp11_1 = tmp1_0*w39;
+                                     const double tmp11_1 = tmp1_0*w39;
-                                     const register double tmp52_1 = B_1_1*w29;
+                                     const double tmp52_1 = B_1_1*w29;
-                                     const register double tmp44_1 = B_1_3*w32;
+                                     const double tmp44_1 = B_1_3*w32;
-                                     const register double tmp25_1 = B_1_1*w39;
+                                     const double tmp25_1 = B_1_1*w39;
-                                     const register double tmp47_1 = tmp2_0*w33;
+                                     const double tmp47_1 = tmp2_0*w33;
-                                     const register double tmp72_1 = B_1_3*w29;
+                                     const double tmp72_1 = B_1_3*w29;
-                                     const register double tmp40_1 = tmp2_0*w28;
+                                     const double tmp40_1 = tmp2_0*w28;
-                                     const register double tmp46_1 = B_1_2*w40;
+                                     const double tmp46_1 = B_1_2*w40;
-                                     const register double tmp36_1 = B_1_2*w42;
+                                     const double tmp36_1 = B_1_2*w42;
-                                     const register double tmp24_1 = B_0_0*w37;
+                                     const double tmp24_1 = B_0_0*w37;
-                                     const register double tmp77_1 = B_0_3*w35;
+                                     const double tmp77_1 = B_0_3*w35;
-                                     const register double tmp68_1 = B_0_3*w37;
+                                     const double tmp68_1 = B_0_3*w37;
-                                     const register double tmp78_1 = B_0_0*w34;
+                                     const double tmp78_1 = B_0_0*w34;
-                                     const register double tmp12_1 = tmp0_0*w36;
+                                     const double tmp12_1 = tmp0_0*w36;
-                                     const register double tmp75_1 = B_1_0*w32;
+                                     const double tmp75_1 = B_1_0*w32;
                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
-Line 2568 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2743 
 void Rectangle::assemblePDESystem(Paso_S
                                      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;
                                  }
                              }
-                         } else { /* constant data */
+                         } else { // constant data
                              for (index_t k=0; k<numEq; k++) {
                                  for (index_t m=0; m<numComp; m++) {
-                                     const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
+                                     const double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
-                                     const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
+                                     const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
-                                     const register double tmp6_1 = B_1*w50;
+                                     const double tmp6_1 = B_1*w50;
-                                     const register double tmp1_1 = B_1*w45;
+                                     const double tmp1_1 = B_1*w45;
-                                     const register double tmp5_1 = B_1*w49;
+                                     const double tmp5_1 = B_1*w49;
-                                     const register double tmp4_1 = B_1*w48;
+                                     const double tmp4_1 = B_1*w48;
-                                     const register double tmp0_1 = B_0*w44;
+                                     const double tmp0_1 = B_0*w44;
-                                     const register double tmp2_1 = B_0*w46;
+                                     const double tmp2_1 = B_0*w46;
-                                     const register double tmp7_1 = B_0*w51;
+                                     const double tmp7_1 = B_0*w51;
-                                     const register double tmp3_1 = B_0*w47;
+                                     const double tmp3_1 = B_0*w47;
                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
-Line 2610 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2785 
 void Rectangle::assemblePDESystem(Paso_S
                          if (C.actsExpanded()) {
                              for (index_t k=0; k<numEq; k++) {
                                  for (index_t m=0; m<numComp; m++) {
-                                     const register double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
+                                     const double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
-                                     const register double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
+                                     const double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
-                                     const register double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
+                                     const double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
-                                     const register double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
+                                     const double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
-                                     const register double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];
+                                     const double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];
-                                     const register double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
+                                     const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
-                                     const register double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];
+                                     const double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];
-                                     const register double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];
+                                     const double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];
-                                     const register double tmp2_0 = C_0_1 + C_0_3;
+                                     const double tmp2_0 = C_0_1 + C_0_3;
-                                     const register double tmp1_0 = C_1_2 + C_1_3;
+                                     const double tmp1_0 = C_1_2 + C_1_3;
-                                     const register double tmp3_0 = C_0_0 + C_0_2;
+                                     const double tmp3_0 = C_0_0 + C_0_2;
-                                     const register double tmp0_0 = C_1_0 + C_1_1;
+                                     const double tmp0_0 = C_1_0 + C_1_1;
-                                     const register double tmp64_1 = C_0_2*w30;
+                                     const double tmp64_1 = C_0_2*w30;
-                                     const register double tmp14_1 = C_0_2*w28;
+                                     const double tmp14_1 = C_0_2*w28;
-                                     const register double tmp19_1 = C_0_3*w31;
+                                     const double tmp19_1 = C_0_3*w31;
-                                     const register double tmp22_1 = C_1_0*w40;
+                                     const double tmp22_1 = C_1_0*w40;
-                                     const register double tmp37_1 = tmp3_0*w35;
+                                     const double tmp37_1 = tmp3_0*w35;
-                                     const register double tmp29_1 = C_0_1*w35;
+                                     const double tmp29_1 = C_0_1*w35;
-                                     const register double tmp73_1 = C_0_2*w37;
+                                     const double tmp73_1 = C_0_2*w37;
-                                     const register double tmp74_1 = C_1_2*w41;
+                                     const double tmp74_1 = C_1_2*w41;
-                                     const register double tmp52_1 = C_1_3*w39;
+                                     const double tmp52_1 = C_1_3*w39;
-                                     const register double tmp25_1 = C_1_1*w39;
+                                     const double tmp25_1 = C_1_1*w39;
-                                     const register double tmp62_1 = C_0_1*w31;
+                                     const double tmp62_1 = C_0_1*w31;
-                                     const register double tmp79_1 = C_1_1*w40;
+                                     const double tmp79_1 = C_1_1*w40;
-                                     const register double tmp43_1 = C_1_1*w36;
+                                     const double tmp43_1 = C_1_1*w36;
-                                     const register double tmp27_1 = C_0_3*w38;
+                                     const double tmp27_1 = C_0_3*w38;
-                                     const register double tmp28_1 = C_1_0*w42;
+                                     const double tmp28_1 = C_1_0*w42;
-                                     const register double tmp63_1 = C_1_1*w42;
+                                     const double tmp63_1 = C_1_1*w42;
-                                     const register double tmp59_1 = C_0_3*w34;
+                                     const double tmp59_1 = C_0_3*w34;
-                                     const register double tmp72_1 = C_1_3*w29;
+                                     const double tmp72_1 = C_1_3*w29;
-                                     const register double tmp40_1 = tmp2_0*w35;
+                                     const double tmp40_1 = tmp2_0*w35;
-                                     const register double tmp13_1 = C_0_3*w30;
+                                     const double tmp13_1 = C_0_3*w30;
-                                     const register double tmp51_1 = C_1_2*w40;
+                                     const double tmp51_1 = C_1_2*w40;
-                                     const register double tmp54_1 = C_1_2*w42;
+                                     const double tmp54_1 = C_1_2*w42;
-                                     const register double tmp12_1 = C_0_0*w31;
+                                     const double tmp12_1 = C_0_0*w31;
-                                     const register double tmp2_1 = tmp1_0*w32;
+                                     const double tmp2_1 = tmp1_0*w32;
-                                     const register double tmp68_1 = C_0_2*w31;
+                                     const double tmp68_1 = C_0_2*w31;
-                                     const register double tmp75_1 = C_1_0*w32;
+                                     const double tmp75_1 = C_1_0*w32;
-                                     const register double tmp49_1 = C_1_1*w41;
+                                     const double tmp49_1 = C_1_1*w41;
-                                     const register double tmp4_1 = C_0_2*w35;
+                                     const double tmp4_1 = C_0_2*w35;
-                                     const register double tmp66_1 = C_0_3*w28;
+                                     const double tmp66_1 = C_0_3*w28;
-                                     const register double tmp56_1 = C_0_1*w37;
+                                     const double tmp56_1 = C_0_1*w37;
-                                     const register double tmp5_1 = C_0_1*w34;
+                                     const double tmp5_1 = C_0_1*w34;
-                                     const register double tmp38_1 = tmp2_0*w34;
+                                     const double tmp38_1 = tmp2_0*w34;
-                                     const register double tmp76_1 = C_0_1*w38;
+                                     const double tmp76_1 = C_0_1*w38;
-                                     const register double tmp21_1 = C_0_1*w28;
+                                     const double tmp21_1 = C_0_1*w28;
-                                     const register double tmp69_1 = C_0_1*w30;
+                                     const double tmp69_1 = C_0_1*w30;
-                                     const register double tmp53_1 = C_1_0*w36;
+                                     const double tmp53_1 = C_1_0*w36;
-                                     const register double tmp42_1 = C_1_2*w39;
+                                     const double tmp42_1 = C_1_2*w39;
-                                     const register double tmp32_1 = tmp1_0*w29;
+                                     const double tmp32_1 = tmp1_0*w29;
-                                     const register double tmp45_1 = C_1_0*w43;
+                                     const double tmp45_1 = C_1_0*w43;
-                                     const register double tmp33_1 = tmp0_0*w32;
+                                     const double tmp33_1 = tmp0_0*w32;
-                                     const register double tmp35_1 = C_1_0*w41;
+                                     const double tmp35_1 = C_1_0*w41;
-                                     const register double tmp26_1 = C_1_2*w36;
+                                     const double tmp26_1 = C_1_2*w36;
-                                     const register double tmp67_1 = C_0_0*w33;
+                                     const double tmp67_1 = C_0_0*w33;
-                                     const register double tmp31_1 = C_0_2*w34;
+                                     const double tmp31_1 = C_0_2*w34;
-                                     const register double tmp20_1 = C_0_0*w30;
+                                     const double tmp20_1 = C_0_0*w30;
-                                     const register double tmp70_1 = C_0_0*w28;
+                                     const double tmp70_1 = C_0_0*w28;
-                                     const register double tmp8_1 = tmp0_0*w39;
+                                     const double tmp8_1 = tmp0_0*w39;
-                                     const register double tmp30_1 = C_1_3*w43;
+                                     const double tmp30_1 = C_1_3*w43;
-                                     const register double tmp0_1 = tmp0_0*w29;
+                                     const double tmp0_1 = tmp0_0*w29;
-                                     const register double tmp17_1 = C_1_3*w41;
+                                     const double tmp17_1 = C_1_3*w41;
-                                     const register double tmp58_1 = C_0_0*w35;
+                                     const double tmp58_1 = C_0_0*w35;
-                                     const register double tmp9_1 = tmp3_0*w33;
+                                     const double tmp9_1 = tmp3_0*w33;
-                                     const register double tmp61_1 = C_1_3*w36;
+                                     const double tmp61_1 = C_1_3*w36;
-                                     const register double tmp41_1 = tmp3_0*w34;
+                                     const double tmp41_1 = tmp3_0*w34;
-                                     const register double tmp50_1 = C_1_3*w32;
+                                     const double tmp50_1 = C_1_3*w32;
-                                     const register double tmp18_1 = C_1_1*w32;
+                                     const double tmp18_1 = C_1_1*w32;
-                                     const register double tmp6_1 = tmp1_0*w36;
+                                     const double tmp6_1 = tmp1_0*w36;
-                                     const register double tmp3_1 = C_0_0*w38;
+                                     const double tmp3_1 = C_0_0*w38;
-                                     const register double tmp34_1 = C_1_1*w29;
+                                     const double tmp34_1 = C_1_1*w29;
-                                     const register double tmp77_1 = C_0_3*w35;
+                                     const double tmp77_1 = C_0_3*w35;
-                                     const register double tmp65_1 = C_1_2*w43;
+                                     const double tmp65_1 = C_1_2*w43;
-                                     const register double tmp71_1 = C_0_3*w33;
+                                     const double tmp71_1 = C_0_3*w33;
-                                     const register double tmp55_1 = C_1_1*w43;
+                                     const double tmp55_1 = C_1_1*w43;
-                                     const register double tmp46_1 = tmp3_0*w28;
+                                     const double tmp46_1 = tmp3_0*w28;
-                                     const register double tmp24_1 = C_0_0*w37;
+                                     const double tmp24_1 = C_0_0*w37;
-                                     const register double tmp10_1 = tmp1_0*w39;
+                                     const double tmp10_1 = tmp1_0*w39;
-                                     const register double tmp48_1 = C_1_0*w29;
+                                     const double tmp48_1 = C_1_0*w29;
-                                     const register double tmp15_1 = C_0_1*w33;
+                                     const double tmp15_1 = C_0_1*w33;
-                                     const register double tmp36_1 = C_1_2*w32;
+                                     const double tmp36_1 = C_1_2*w32;
-                                     const register double tmp60_1 = C_1_0*w39;
+                                     const double tmp60_1 = C_1_0*w39;
-                                     const register double tmp47_1 = tmp2_0*w33;
+                                     const double tmp47_1 = tmp2_0*w33;
-                                     const register double tmp16_1 = C_1_2*w29;
+                                     const double tmp16_1 = C_1_2*w29;
-                                     const register double tmp1_1 = C_0_3*w37;
+                                     const double tmp1_1 = C_0_3*w37;
-                                     const register double tmp7_1 = tmp2_0*w28;
+                                     const double tmp7_1 = tmp2_0*w28;
-                                     const register double tmp39_1 = C_1_3*w40;
+                                     const double tmp39_1 = C_1_3*w40;
-                                     const register double tmp44_1 = C_1_3*w42;
+                                     const double tmp44_1 = C_1_3*w42;
-                                     const register double tmp57_1 = C_0_2*w38;
+                                     const double tmp57_1 = C_0_2*w38;
-                                     const register double tmp78_1 = C_0_0*w34;
+                                     const double tmp78_1 = C_0_0*w34;
-                                     const register double tmp11_1 = tmp0_0*w36;
+                                     const double tmp11_1 = tmp0_0*w36;
-                                     const register double tmp23_1 = C_0_2*w33;
+                                     const double tmp23_1 = C_0_2*w33;
                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
-Line 2720 
 void Rectangle::assemblePDESystem(Paso_S
+Line 2895 
 void Rectangle::assemblePDESystem(Paso_S
                                      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;
                                  }
                              }
-                         } else { /* constant data */
+                         } else { // constant data
                              for (index_t k=0; k<numEq; k++) {
                                  for (index_t m=0; m<numComp; m++) {
-                                     const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
+                                     const double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
-                                     const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
+                                     const double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
-                                     const register double tmp1_1 = C_1*w45;
+                                     const double tmp1_1 = C_1*w45;
-                                     const register double tmp3_1 = C_0*w51;
+                                     const double tmp3_1 = C_0*w51;
-                                     const register double tmp4_1 = C_0*w44;
+                                     const double tmp4_1 = C_0*w44;