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

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

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

revision 3760 by caltinay, Mon Jan 9 05:21:18 2012 UTC revision 3770 by gross, Wed Jan 18 01:40:15 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2011 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 49  Rectangle::Rectangle(int n0, int n1, dou Line 49  Rectangle::Rectangle(int n0, int n1, dou
49      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
50          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
51    
52      // local number of elements (including overlap)      // local number of elements (with and without overlap)
53      m_NE0 = (m_NX>1 ? (n0+1)/m_NX : n0);      m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
54      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
55          m_NE0++;          m_NE0++;
56      m_NE1 = (m_NY>1 ? (n1+1)/m_NY : n1);      else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)
57            m_ownNE0--;
58    
59        m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
60      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)
61          m_NE1++;          m_NE1++;
62        else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)
63            m_ownNE1--;
64    
65      // local number of nodes      // local number of nodes
66      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
# Line 106  void Rectangle::dump(const string& fileN Line 111  void Rectangle::dump(const string& fileN
111      const int NUM_SILO_FILES = 1;      const int NUM_SILO_FILES = 1;
112      const char* blockDirFmt = "/block%04d";      const char* blockDirFmt = "/block%04d";
113      int driver=DB_HDF5;          int driver=DB_HDF5;    
     string siloPath;  
114      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
115    
116  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 126  void Rectangle::dump(const string& fileN Line 130  void Rectangle::dump(const string& fileN
130                          PMPIO_DefaultClose, (void*)&driver);                          PMPIO_DefaultClose, (void*)&driver);
131          }          }
132          if (baton) {          if (baton) {
133              char str[64];              char siloPath[64];
134              snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));              snprintf(siloPath, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
135              siloPath = str;              dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath);
             dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());  
136          }          }
137  #endif  #endif
138      } else {      } else {
# Line 141  void Rectangle::dump(const string& fileN Line 144  void Rectangle::dump(const string& fileN
144              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,              dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
145                      getDescription().c_str(), driver);                      getDescription().c_str(), driver);
146          }          }
147            char siloPath[64];
148            snprintf(siloPath, 64, blockDirFmt, 0);
149            DBMkDir(dbfile, siloPath);
150            DBSetDir(dbfile, siloPath);
151      }      }
152    
153      if (!dbfile)      if (!dbfile)
# Line 241  const int* Rectangle::borrowSampleRefere Line 248  const int* Rectangle::borrowSampleRefere
248  {  {
249      switch (fsType) {      switch (fsType) {
250          case Nodes:          case Nodes:
251          case ReducedNodes: //FIXME: reduced          case ReducedNodes: // FIXME: reduced
252              return &m_nodeId[0];              return &m_nodeId[0];
253          case DegreesOfFreedom:          case DegreesOfFreedom:
254          case ReducedDegreesOfFreedom: //FIXME: reduced          case ReducedDegreesOfFreedom: // FIXME: reduced
255              return &m_dofId[0];              return &m_dofId[0];
256          case Elements:          case Elements:
257          case ReducedElements:          case ReducedElements:
# Line 269  bool Rectangle::ownSample(int fsType, in Line 276  bool Rectangle::ownSample(int fsType, in
276    
277      switch (fsType) {      switch (fsType) {
278          case Nodes:          case Nodes:
279          case ReducedNodes: //FIXME: reduced          case ReducedNodes: // FIXME: reduced
280              return (m_dofMap[id] < getNumDOF());              return (m_dofMap[id] < getNumDOF());
281          case DegreesOfFreedom:          case DegreesOfFreedom:
282          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
# Line 281  bool Rectangle::ownSample(int fsType, in Line 288  bool Rectangle::ownSample(int fsType, in
288          case FaceElements:          case FaceElements:
289          case ReducedFaceElements:          case ReducedFaceElements:
290              {              {
291                  // check ownership of face element's first node                  // determine which face the sample belongs to before
292                    // checking ownership of corresponding element's first node
293                  const IndexVector faces = getNumFacesPerBoundary();                  const IndexVector faces = getNumFacesPerBoundary();
294                  dim_t n=0;                  dim_t n=0;
295                  for (size_t i=0; i<faces.size(); i++) {                  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
297                      if (id<n) {                      if (id<n) {
298                          index_t k;                          index_t k;
299                          if (i==1)                          if (i==1)
300                              k=m_N0-1;                              k=m_N0-2;
301                          else if (i==3)                          else if (i==3)
302                              k=m_N0*(m_N1-1);                              k=m_N0*(m_N1-2);
303                          else                          else
304                              k=0;                              k=0;
305                          // determine whether to move right or up                          // determine whether to move right or up
# Line 311  bool Rectangle::ownSample(int fsType, in Line 319  bool Rectangle::ownSample(int fsType, in
319      throw RipleyException(msg.str());      throw RipleyException(msg.str());
320  }  }
321    
 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());  
     }  
 }  
   
322  void Rectangle::setToNormal(escript::Data& out) const  void Rectangle::setToNormal(escript::Data& out) const
323  {  {
324      if (out.getFunctionSpace().getTypeCode() == FaceElements) {      if (out.getFunctionSpace().getTypeCode() == FaceElements) {
# Line 906  IndexVector Rectangle::getNumFacesPerBou Line 545  IndexVector Rectangle::getNumFacesPerBou
545      return ret;      return ret;
546  }  }
547    
548    IndexVector Rectangle::getNumSubdivisionsPerDim() const
549    {
550        IndexVector ret;
551        ret.push_back(m_NX);
552        ret.push_back(m_NY);
553        return ret;
554    }
555    
556  pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
557  {  {
558      if (dim==0) {      if (dim==0) {
# Line 956  void Rectangle::assembleCoordinates(escr Line 603  void Rectangle::assembleCoordinates(escr
603  }  }
604    
605  //protected  //protected
606    void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
607    {
608        const dim_t numComp = in.getDataPointSize();
609        const double h0 = m_l0/m_gNE0;
610        const double h1 = m_l1/m_gNE1;
611        const double cx0 = -1./h0;
612        const double cx1 = -.78867513459481288225/h0;
613        const double cx2 = -.5/h0;
614        const double cx3 = -.21132486540518711775/h0;
615        const double cx4 = .21132486540518711775/h0;
616        const double cx5 = .5/h0;
617        const double cx6 = .78867513459481288225/h0;
618        const double cx7 = 1./h0;
619        const double cy0 = -1./h1;
620        const double cy1 = -.78867513459481288225/h1;
621        const double cy2 = -.5/h1;
622        const double cy3 = -.21132486540518711775/h1;
623        const double cy4 = .21132486540518711775/h1;
624        const double cy5 = .5/h1;
625        const double cy6 = .78867513459481288225/h1;
626        const double cy7 = 1./h1;
627    
628        if (out.getFunctionSpace().getTypeCode() == Elements) {
629            out.requireWrite();
630    #pragma omp parallel for
631            for (index_t k1=0; k1 < m_NE1; ++k1) {
632                for (index_t k0=0; k0 < m_NE0; ++k0) {
633                    const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
634                    const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
635                    const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
636                    const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
637                    double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
638                    for (index_t i=0; i < numComp; ++i) {
639                        o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
640                        o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
641                        o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
642                        o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
643                        o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
644                        o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
645                        o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
646                        o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
647                    } /* end of component loop i */
648                } /* end of k0 loop */
649            } /* end of k1 loop */
650        } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
651            out.requireWrite();
652    #pragma omp parallel for
653            for (index_t k1=0; k1 < m_NE1; ++k1) {
654                for (index_t k0=0; k0 < m_NE0; ++k0) {
655                    const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
656                    const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
657                    const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
658                    const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
659                    double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
660                    for (index_t i=0; i < numComp; ++i) {
661                        o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
662                        o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
663                    } /* end of component loop i */
664                } /* end of k0 loop */
665            } /* end of k1 loop */
666        } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
667            out.requireWrite();
668    #pragma omp parallel
669            {
670                if (m_faceOffset[0] > -1) {
671    #pragma omp for nowait
672                    for (index_t k1=0; k1 < m_NE1; ++k1) {
673                        const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
674                        const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
675                        const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
676                        const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
677                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
678                        for (index_t i=0; i < numComp; ++i) {
679                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
680                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
681                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
682                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
683                        } /* end of component loop i */
684                    } /* end of k1 loop */
685                } /* end of face 0 */
686                if (m_faceOffset[1] > -1) {
687    #pragma omp for nowait
688                    for (index_t k1=0; k1 < m_NE1; ++k1) {
689                        const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
690                        const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
691                        const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
692                        const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
693                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
694                        for (index_t i=0; i < numComp; ++i) {
695                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
696                            o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
697                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
698                            o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
699                        } /* end of component loop i */
700                    } /* end of k1 loop */
701                } /* end of face 1 */
702                if (m_faceOffset[2] > -1) {
703    #pragma omp for nowait
704                    for (index_t k0=0; k0 < m_NE0; ++k0) {
705                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
706                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
707                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
708                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
709                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
710                        for (index_t i=0; i < numComp; ++i) {
711                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
712                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
713                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
714                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
715                        } /* end of component loop i */
716                    } /* end of k0 loop */
717                } /* end of face 2 */
718                if (m_faceOffset[3] > -1) {
719    #pragma omp for nowait
720                    for (index_t k0=0; k0 < m_NE0; ++k0) {
721                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
722                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
723                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
724                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
725                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
726                        for (index_t i=0; i < numComp; ++i) {
727                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
728                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
729                            o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
730                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
731                        } /* end of component loop i */
732                    } /* end of k0 loop */
733                } /* end of face 3 */
734            } // end of parallel section
735        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
736            out.requireWrite();
737    #pragma omp parallel
738            {
739                if (m_faceOffset[0] > -1) {
740    #pragma omp for nowait
741                    for (index_t k1=0; k1 < m_NE1; ++k1) {
742                        const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
743                        const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
744                        const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
745                        const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
746                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
747                        for (index_t i=0; i < numComp; ++i) {
748                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
749                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
750                        } /* end of component loop i */
751                    } /* end of k1 loop */
752                } /* end of face 0 */
753                if (m_faceOffset[1] > -1) {
754    #pragma omp for nowait
755                    for (index_t k1=0; k1 < m_NE1; ++k1) {
756                        const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
757                        const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
758                        const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
759                        const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
760                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
761                        for (index_t i=0; i < numComp; ++i) {
762                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
763                            o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
764                        } /* end of component loop i */
765                    } /* end of k1 loop */
766                } /* end of face 1 */
767                if (m_faceOffset[2] > -1) {
768    #pragma omp for nowait
769                    for (index_t k0=0; k0 < m_NE0; ++k0) {
770                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
771                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
772                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
773                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
774                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
775                        for (index_t i=0; i < numComp; ++i) {
776                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
777                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
778                        } /* end of component loop i */
779                    } /* end of k0 loop */
780                } /* end of face 2 */
781                if (m_faceOffset[3] > -1) {
782    #pragma omp for nowait
783                    for (index_t k0=0; k0 < m_NE0; ++k0) {
784                        const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
785                        const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
786                        const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
787                        const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
788                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
789                        for (index_t i=0; i < numComp; ++i) {
790                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
791                            o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
792                        } /* end of component loop i */
793                    } /* end of k0 loop */
794                } /* end of face 3 */
795            } // end of parallel section
796        }
797    }
798    
799    //protected
800    void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
801    {
802        const dim_t numComp = arg.getDataPointSize();
803        const double h0 = m_l0/m_gNE0;
804        const double h1 = m_l1/m_gNE1;
805        const index_t left = (m_offset0==0 ? 0 : 1);
806        const index_t bottom = (m_offset1==0 ? 0 : 1);
807        if (arg.getFunctionSpace().getTypeCode() == Elements) {
808            const double w = h0*h1/4.;
809    #pragma omp parallel
810            {
811                vector<double> int_local(numComp, 0);
812    #pragma omp for nowait
813                for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
814                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
815                        const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
816                        for (index_t i=0; i < numComp; ++i) {
817                            const double f0 = f[INDEX2(i,0,numComp)];
818                            const double f1 = f[INDEX2(i,1,numComp)];
819                            const double f2 = f[INDEX2(i,2,numComp)];
820                            const double f3 = f[INDEX2(i,3,numComp)];
821                            int_local[i]+=(f0+f1+f2+f3)*w;
822                        }  /* end of component loop i */
823                    } /* end of k0 loop */
824                } /* end of k1 loop */
825    
826    #pragma omp critical
827                for (index_t i=0; i<numComp; i++)
828                    integrals[i]+=int_local[i];
829            } // end of parallel section
830        } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
831            const double w = h0*h1;
832    #pragma omp parallel
833            {
834                vector<double> int_local(numComp, 0);
835    #pragma omp for nowait
836                for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
837                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
838                        const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
839                        for (index_t i=0; i < numComp; ++i) {
840                            int_local[i]+=f[i]*w;
841                        }  /* end of component loop i */
842                    } /* end of k0 loop */
843                } /* end of k1 loop */
844    
845    #pragma omp critical
846                for (index_t i=0; i<numComp; i++)
847                    integrals[i]+=int_local[i];
848            } // end of parallel section
849        } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
850            const double w0 = h0/2.;
851            const double w1 = h1/2.;
852    #pragma omp parallel
853            {
854                vector<double> int_local(numComp, 0);
855                if (m_faceOffset[0] > -1) {
856    #pragma omp for nowait
857                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
858                        const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
859                        for (index_t i=0; i < numComp; ++i) {
860                            const double f0 = f[INDEX2(i,0,numComp)];
861                            const double f1 = f[INDEX2(i,1,numComp)];
862                            int_local[i]+=(f0+f1)*w1;
863                        }  /* end of component loop i */
864                    } /* end of k1 loop */
865                }
866    
867                if (m_faceOffset[1] > -1) {
868    #pragma omp for nowait
869                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
870                        const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
871                        for (index_t i=0; i < numComp; ++i) {
872                            const double f0 = f[INDEX2(i,0,numComp)];
873                            const double f1 = f[INDEX2(i,1,numComp)];
874                            int_local[i]+=(f0+f1)*w1;
875                        }  /* end of component loop i */
876                    } /* end of k1 loop */
877                }
878    
879                if (m_faceOffset[2] > -1) {
880    #pragma omp for nowait
881                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
882                        const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
883                        for (index_t i=0; i < numComp; ++i) {
884                            const double f0 = f[INDEX2(i,0,numComp)];
885                            const double f1 = f[INDEX2(i,1,numComp)];
886                            int_local[i]+=(f0+f1)*w0;
887                        }  /* end of component loop i */
888                    } /* end of k0 loop */
889                }
890    
891                if (m_faceOffset[3] > -1) {
892    #pragma omp for nowait
893                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
894                        const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
895                        for (index_t i=0; i < numComp; ++i) {
896                            const double f0 = f[INDEX2(i,0,numComp)];
897                            const double f1 = f[INDEX2(i,1,numComp)];
898                            int_local[i]+=(f0+f1)*w0;
899                        }  /* end of component loop i */
900                    } /* end of k0 loop */
901                }
902    
903    #pragma omp critical
904                for (index_t i=0; i<numComp; i++)
905                    integrals[i]+=int_local[i];
906            } // end of parallel section
907        } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
908    #pragma omp parallel
909            {
910                vector<double> int_local(numComp, 0);
911                if (m_faceOffset[0] > -1) {
912    #pragma omp for nowait
913                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
914                        const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
915                        for (index_t i=0; i < numComp; ++i) {
916                            int_local[i]+=f[i]*h1;
917                        }  /* end of component loop i */
918                    } /* end of k1 loop */
919                }
920    
921                if (m_faceOffset[1] > -1) {
922    #pragma omp for nowait
923                    for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
924                        const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
925                        for (index_t i=0; i < numComp; ++i) {
926                            int_local[i]+=f[i]*h1;
927                        }  /* end of component loop i */
928                    } /* end of k1 loop */
929                }
930    
931                if (m_faceOffset[2] > -1) {
932    #pragma omp for nowait
933                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
934                        const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
935                        for (index_t i=0; i < numComp; ++i) {
936                            int_local[i]+=f[i]*h0;
937                        }  /* end of component loop i */
938                    } /* end of k0 loop */
939                }
940    
941                if (m_faceOffset[3] > -1) {
942    #pragma omp for nowait
943                    for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
944                        const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
945                        for (index_t i=0; i < numComp; ++i) {
946                            int_local[i]+=f[i]*h0;
947                        }  /* end of component loop i */
948                    } /* end of k0 loop */
949                }
950    
951    #pragma omp critical
952                for (index_t i=0; i<numComp; i++)
953                    integrals[i]+=int_local[i];
954            } // end of parallel section
955        }
956    }
957    
958    //protected
959  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
960  {  {
961      const dim_t nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
# Line 1110  void Rectangle::createPattern() Line 1110  void Rectangle::createPattern()
1110      // The rest is assigned in the loop further down      // The rest is assigned in the loop further down
1111      m_dofMap.assign(getNumNodes(), 0);      m_dofMap.assign(getNumNodes(), 0);
1112  #pragma omp parallel for  #pragma omp parallel for
1113      for (index_t i=bottom; i<m_N1; i++) {      for (index_t i=bottom; i<bottom+nDOF1; i++) {
1114          for (index_t j=left; j<m_N0; j++) {          for (index_t j=left; j<left+nDOF0; j++) {
1115              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
1116          }          }
1117      }      }
# Line 1279  void Rectangle::interpolateNodesOnElemen Line 1279  void Rectangle::interpolateNodesOnElemen
1279      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1280      if (reduced) {      if (reduced) {
1281          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */  
1282          const double c0 = .25;          const double c0 = .25;
1283  #pragma omp parallel for  #pragma omp parallel for
1284          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1285              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
1286                  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));
1287                  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));
1288                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1289                  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));
1290                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1291                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1292                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1293                  } /* end of component loop i */                  } /* end of component loop i */
1294              } /* end of k0 loop */              } /* end of k0 loop */
1295          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */  
1296      } else {      } else {
1297          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */  
1298          const double c0 = .16666666666666666667;          const double c0 = .16666666666666666667;
1299          const double c1 = .044658198738520451079;          const double c1 = .044658198738520451079;
1300          const double c2 = .62200846792814621559;          const double c2 = .62200846792814621559;
1301  #pragma omp parallel for  #pragma omp parallel for
1302          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1303              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
1304                  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));
1305                  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));
1306                  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));
1307                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1308                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1309                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1310                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);
# Line 1317  void Rectangle::interpolateNodesOnElemen Line 1314  void Rectangle::interpolateNodesOnElemen
1314                  } /* end of component loop i */                  } /* end of component loop i */
1315              } /* end of k0 loop */              } /* end of k0 loop */
1316          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */  
1317      }      }
1318  }  }
1319    
# Line 1331  void Rectangle::interpolateNodesOnFaces( Line 1327  void Rectangle::interpolateNodesOnFaces(
1327          const double c0 = .5;          const double c0 = .5;
1328  #pragma omp parallel  #pragma omp parallel
1329          {          {
             /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */  
1330              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1331  #pragma omp for nowait  #pragma omp for nowait
1332                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1333                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1334                      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));
1335                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1336                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1337                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
# Line 1346  void Rectangle::interpolateNodesOnFaces( Line 1341  void Rectangle::interpolateNodesOnFaces(
1341              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1342  #pragma omp for nowait  #pragma omp for nowait
1343                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1344                      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));
1345                      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));
1346                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1347                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1348                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
# Line 1357  void Rectangle::interpolateNodesOnFaces( Line 1352  void Rectangle::interpolateNodesOnFaces(
1352              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1353  #pragma omp for nowait  #pragma omp for nowait
1354                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1355                      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));
1356                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1357                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1358                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1359                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
# Line 1368  void Rectangle::interpolateNodesOnFaces( Line 1363  void Rectangle::interpolateNodesOnFaces(
1363              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1364  #pragma omp for nowait  #pragma omp for nowait
1365                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1366                      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));
1367                      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));
1368                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1369                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1370                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1371                      } /* end of component loop i */                      } /* end of component loop i */
1372                  } /* end of k0 loop */                  } /* end of k0 loop */
1373              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */  
1374          } // end of parallel section          } // end of parallel section
1375      } else {      } else {
1376          out.requireWrite();          out.requireWrite();
# Line 1384  void Rectangle::interpolateNodesOnFaces( Line 1378  void Rectangle::interpolateNodesOnFaces(
1378          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1379  #pragma omp parallel  #pragma omp parallel
1380          {          {
             /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */  
1381              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1382  #pragma omp for nowait  #pragma omp for nowait
1383                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1384                      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));
1385                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1386                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1387                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1388                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;
# Line 1400  void Rectangle::interpolateNodesOnFaces( Line 1393  void Rectangle::interpolateNodesOnFaces(
1393              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1394  #pragma omp for nowait  #pragma omp for nowait
1395                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1396                      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));
1397                      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));
1398                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1399                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1400                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;
# Line 1412  void Rectangle::interpolateNodesOnFaces( Line 1405  void Rectangle::interpolateNodesOnFaces(
1405              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1406  #pragma omp for nowait  #pragma omp for nowait
1407                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1408                      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));
1409                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1410                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1411                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1412                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;
# Line 1424  void Rectangle::interpolateNodesOnFaces( Line 1417  void Rectangle::interpolateNodesOnFaces(
1417              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1418  #pragma omp for nowait  #pragma omp for nowait
1419                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1420                      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));
1421                      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));
1422                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1423                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1424                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;
# Line 1433  void Rectangle::interpolateNodesOnFaces( Line 1426  void Rectangle::interpolateNodesOnFaces(
1426                      } /* end of component loop i */                      } /* end of component loop i */
1427                  } /* end of k0 loop */                  } /* end of k0 loop */
1428              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */  
1429          } // end of parallel section          } // end of parallel section
1430      }      }
1431  }  }
# Line 1442  void Rectangle::interpolateNodesOnFaces( Line 1434  void Rectangle::interpolateNodesOnFaces(
1434  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1435          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1436          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
1437          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  
1438  {  {
1439      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1440      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
     /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */  
1441      const double w0 = -0.1555021169820365539*h1/h0;      const double w0 = -0.1555021169820365539*h1/h0;
1442      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
1443        const double w2 = -0.15550211698203655390;
1444        const double w3 = 0.041666666666666666667*h0/h1;
1445        const double w4 = 0.15550211698203655390;
1446        const double w5 = -0.041666666666666666667;
1447        const double w6 = -0.01116454968463011277*h1/h0;
1448        const double w7 = 0.011164549684630112770;
1449        const double w8 = -0.011164549684630112770;
1450        const double w9 = -0.041666666666666666667*h1/h0;
1451      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.041666666666666666667*h0/h1;
1452      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = 0.1555021169820365539*h1/h0;
1453      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = 0.1555021169820365539*h0/h1;
# Line 1459  void Rectangle::assemblePDESingle(Paso_S Line 1457  void Rectangle::assemblePDESingle(Paso_S
1457      const double w16 = -0.01116454968463011277*h0/h1;      const double w16 = -0.01116454968463011277*h0/h1;
1458      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*h0/h1;
1459      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*h1/h0;
1460      const double w19 = 0.25000000000000000000;      const double w19 = 0.25;
1461      const double w2 = -0.15550211698203655390;      const double w20 = -0.25;
     const double w20 = -0.25000000000000000000;  
1462      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*h0/h1;
1463      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*h1/h0;
1464      const double w23 = -0.16666666666666666667*h0/h1;      const double w23 = -0.16666666666666666667*h0/h1;
# Line 1471  void Rectangle::assemblePDESingle(Paso_S Line 1468  void Rectangle::assemblePDESingle(Paso_S
1468      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = -0.33333333333333333333*h0/h1;
1469      const double w28 = -0.032861463941450536761*h1;      const double w28 = -0.032861463941450536761*h1;
1470      const double w29 = -0.032861463941450536761*h0;      const double w29 = -0.032861463941450536761*h0;
     const double w3 = 0.041666666666666666667*h0/h1;  
1471      const double w30 = -0.12264065304058601714*h1;      const double w30 = -0.12264065304058601714*h1;
1472      const double w31 = -0.0023593469594139828636*h1;      const double w31 = -0.0023593469594139828636*h1;
1473      const double w32 = -0.008805202725216129906*h0;      const double w32 = -0.008805202725216129906*h0;
# Line 1482  void Rectangle::assemblePDESingle(Paso_S Line 1478  void Rectangle::assemblePDESingle(Paso_S
1478      const double w37 = 0.0023593469594139828636*h1;      const double w37 = 0.0023593469594139828636*h1;
1479      const double w38 = 0.12264065304058601714*h1;      const double w38 = 0.12264065304058601714*h1;
1480      const double w39 = 0.032861463941450536761*h0;      const double w39 = 0.032861463941450536761*h0;
     const double w4 = 0.15550211698203655390;  
1481      const double w40 = -0.12264065304058601714*h0;      const double w40 = -0.12264065304058601714*h0;
1482      const double w41 = -0.0023593469594139828636*h0;      const double w41 = -0.0023593469594139828636*h0;
1483      const double w42 = 0.0023593469594139828636*h0;      const double w42 = 0.0023593469594139828636*h0;
# Line 1493  void Rectangle::assemblePDESingle(Paso_S Line 1488  void Rectangle::assemblePDESingle(Paso_S
1488      const double w47 = 0.16666666666666666667*h1;      const double w47 = 0.16666666666666666667*h1;
1489      const double w48 = 0.083333333333333333333*h0;      const double w48 = 0.083333333333333333333*h0;
1490      const double w49 = -0.16666666666666666667*h0;      const double w49 = -0.16666666666666666667*h0;
     const double w5 = -0.041666666666666666667;  
1491      const double w50 = 0.16666666666666666667*h0;      const double w50 = 0.16666666666666666667*h0;
1492      const double w51 = -0.083333333333333333333*h1;      const double w51 = -0.083333333333333333333*h1;
1493      const double w52 = 0.025917019497006092316*h0*h1;      const double w52 = 0.025917019497006092316*h0*h1;
# Line 1504  void Rectangle::assemblePDESingle(Paso_S Line 1498  void Rectangle::assemblePDESingle(Paso_S
1498      const double w57 = 0.055555555555555555556*h0*h1;      const double w57 = 0.055555555555555555556*h0*h1;
1499      const double w58 = 0.027777777777777777778*h0*h1;      const double w58 = 0.027777777777777777778*h0*h1;
1500      const double w59 = 0.11111111111111111111*h0*h1;      const double w59 = 0.11111111111111111111*h0*h1;
     const double w6 = -0.01116454968463011277*h1/h0;  
1501      const double w60 = -0.19716878364870322056*h1;      const double w60 = -0.19716878364870322056*h1;
1502      const double w61 = -0.19716878364870322056*h0;      const double w61 = -0.19716878364870322056*h0;
1503      const double w62 = -0.052831216351296779436*h0;      const double w62 = -0.052831216351296779436*h0;
# Line 1515  void Rectangle::assemblePDESingle(Paso_S Line 1508  void Rectangle::assemblePDESingle(Paso_S
1508      const double w67 = 0.052831216351296779436*h0;      const double w67 = 0.052831216351296779436*h0;
1509      const double w68 = -0.5*h1;      const double w68 = -0.5*h1;
1510      const double w69 = -0.5*h0;      const double w69 = -0.5*h0;
     const double w7 = 0.011164549684630112770;  
1511      const double w70 = 0.5*h1;      const double w70 = 0.5*h1;
1512      const double w71 = 0.5*h0;      const double w71 = 0.5*h0;
1513      const double w72 = 0.1555021169820365539*h0*h1;      const double w72 = 0.1555021169820365539*h0*h1;
1514      const double w73 = 0.041666666666666666667*h0*h1;      const double w73 = 0.041666666666666666667*h0*h1;
1515      const double w74 = 0.01116454968463011277*h0*h1;      const double w74 = 0.01116454968463011277*h0*h1;
1516      const double w75 = 0.25*h0*h1;      const double w75 = 0.25*h0*h1;
     const double w8 = -0.011164549684630112770;  
     const double w9 = -0.041666666666666666667*h1/h0;  
     /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */  
1517    
1518      rhs.requireWrite();      rhs.requireWrite();
1519  #pragma omp parallel  #pragma omp parallel
# Line 1538  void Rectangle::assemblePDESingle(Paso_S Line 1527  void Rectangle::assemblePDESingle(Paso_S
1527                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1528                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1529                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
                     /*** GENERATOR SNIP_PDE_SINGLE TOP */  
1530                      ///////////////                      ///////////////
1531                      // process A //                      // process A //
1532                      ///////////////                      ///////////////
# Line 1546  void Rectangle::assemblePDESingle(Paso_S Line 1534  void Rectangle::assemblePDESingle(Paso_S
1534                          add_EM_S=true;                          add_EM_S=true;
1535                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1536                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
1537                              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)];
1538                              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)];
1539                              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)];
1540                              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)];
1541                              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)];
1542                              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)];
1543                              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)];
1544                              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)];
1545                              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)];
1546                              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)];
1547                              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)];
1548                              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)];
1549                              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)];
1550                              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)];
1551                              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)];
1552                              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)];
1553                              const register double tmp4_0 = A_10_1 + A_10_2;                              const double tmp0_0 = A_01_0 + A_01_3;
1554                              const register double tmp12_0 = A_11_0 + A_11_2;                              const double tmp1_0 = A_00_0 + A_00_1;
1555                              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;
1556                              const register double tmp10_0 = A_01_3 + A_10_3;                              const double tmp3_0 = A_00_2 + A_00_3;
1557                              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;
1558                              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;
1559                              const register double tmp13_0 = A_01_2 + A_10_1;                              const double tmp6_0 = A_01_3 + A_10_0;
1560                              const register double tmp3_0 = A_00_2 + A_00_3;                              const double tmp7_0 = A_01_0 + A_10_3;
1561                              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;
1562                              const register double tmp18_0 = A_01_1 + A_10_1;                              const double tmp9_0 = A_01_0 + A_10_0;
1563                              const register double tmp1_0 = A_00_0 + A_00_1;                              const double tmp12_0 = A_11_0 + A_11_2;
1564                              const register double tmp15_0 = A_01_1 + A_10_2;                              const double tmp10_0 = A_01_3 + A_10_3;
1565                              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;
1566                              const register double tmp16_0 = A_10_0 + A_10_3;                              const double tmp13_0 = A_01_2 + A_10_1;
1567                              const register double tmp6_0 = A_01_3 + A_10_0;                              const double tmp11_0 = A_11_1 + A_11_3;
1568                              const register double tmp17_0 = A_01_1 + A_01_2;                              const double tmp18_0 = A_01_1 + A_10_1;
1569                              const register double tmp9_0 = A_01_0 + A_10_0;                              const double tmp15_0 = A_01_1 + A_10_2;
1570                              const register double tmp7_0 = A_01_0 + A_10_3;                              const double tmp16_0 = A_10_0 + A_10_3;
1571                              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;
1572                              const register double tmp19_0 = A_01_2 + A_10_2;                              const double tmp19_0 = A_01_2 + A_10_2;
1573                              const register double tmp14_1 = A_10_0*w8;                              const double tmp0_1 = A_10_3*w8;
1574                              const register double tmp23_1 = tmp3_0*w14;                              const double tmp1_1 = tmp0_0*w1;
1575                              const register double tmp35_1 = A_01_0*w8;                              const double tmp2_1 = A_01_1*w4;
1576                              const register double tmp54_1 = tmp13_0*w8;                              const double tmp3_1 = tmp1_0*w0;
1577                              const register double tmp20_1 = tmp9_0*w4;                              const double tmp4_1 = A_01_2*w7;
1578                              const register double tmp25_1 = tmp12_0*w12;                              const double tmp5_1 = tmp2_0*w3;
1579                              const register double tmp2_1 = A_01_1*w4;                              const double tmp6_1 = tmp3_0*w6;
1580                              const register double tmp44_1 = tmp7_0*w7;                              const double tmp7_1 = A_10_0*w2;
1581                              const register double tmp26_1 = tmp10_0*w4;                              const double tmp8_1 = tmp4_0*w5;
1582                              const register double tmp52_1 = tmp18_0*w8;                              const double tmp9_1 = tmp2_0*w10;
1583                              const register double tmp48_1 = A_10_1*w7;                              const double tmp14_1 = A_10_0*w8;
1584                              const register double tmp46_1 = A_01_3*w8;                              const double tmp23_1 = tmp3_0*w14;
1585                              const register double tmp50_1 = A_01_0*w2;                              const double tmp35_1 = A_01_0*w8;
1586                              const register double tmp8_1 = tmp4_0*w5;                              const double tmp54_1 = tmp13_0*w8;
1587                              const register double tmp56_1 = tmp19_0*w8;                              const double tmp20_1 = tmp9_0*w4;
1588                              const register double tmp9_1 = tmp2_0*w10;                              const double tmp25_1 = tmp12_0*w12;
1589                              const register double tmp19_1 = A_10_3*w2;                              const double tmp44_1 = tmp7_0*w7;
1590                              const register double tmp47_1 = A_10_2*w4;                              const double tmp26_1 = tmp10_0*w4;
1591                              const register double tmp16_1 = tmp3_0*w0;                              const double tmp52_1 = tmp18_0*w8;
1592                              const register double tmp18_1 = tmp1_0*w6;                              const double tmp48_1 = A_10_1*w7;
1593                              const register double tmp31_1 = tmp11_0*w12;                              const double tmp46_1 = A_01_3*w8;
1594                              const register double tmp55_1 = tmp15_0*w2;                              const double tmp50_1 = A_01_0*w2;
1595                              const register double tmp39_1 = A_10_2*w7;                              const double tmp56_1 = tmp19_0*w8;
1596                              const register double tmp11_1 = tmp6_0*w7;                              const double tmp19_1 = A_10_3*w2;
1597                              const register double tmp40_1 = tmp11_0*w17;                              const double tmp47_1 = A_10_2*w4;
1598                              const register double tmp34_1 = tmp15_0*w8;                              const double tmp16_1 = tmp3_0*w0;
1599                              const register double tmp33_1 = tmp14_0*w5;                              const double tmp18_1 = tmp1_0*w6;
1600                              const register double tmp24_1 = tmp11_0*w13;                              const double tmp31_1 = tmp11_0*w12;
1601                              const register double tmp3_1 = tmp1_0*w0;                              const double tmp55_1 = tmp15_0*w2;
1602                              const register double tmp5_1 = tmp2_0*w3;                              const double tmp39_1 = A_10_2*w7;
1603                              const register double tmp43_1 = tmp17_0*w5;                              const double tmp11_1 = tmp6_0*w7;
1604                              const register double tmp15_1 = A_01_2*w4;                              const double tmp40_1 = tmp11_0*w17;
1605                              const register double tmp53_1 = tmp19_0*w2;                              const double tmp34_1 = tmp15_0*w8;
1606                              const register double tmp27_1 = tmp3_0*w11;                              const double tmp33_1 = tmp14_0*w5;
1607                              const register double tmp32_1 = tmp13_0*w2;                              const double tmp24_1 = tmp11_0*w13;
1608                              const register double tmp10_1 = tmp5_0*w9;                              const double tmp43_1 = tmp17_0*w5;
1609                              const register double tmp37_1 = A_10_1*w4;                              const double tmp15_1 = A_01_2*w4;
1610                              const register double tmp38_1 = tmp5_0*w15;                              const double tmp53_1 = tmp19_0*w2;
1611                              const register double tmp17_1 = A_01_1*w7;                              const double tmp27_1 = tmp3_0*w11;
1612                              const register double tmp12_1 = tmp7_0*w4;                              const double tmp32_1 = tmp13_0*w2;
1613                              const register double tmp22_1 = tmp10_0*w7;                              const double tmp10_1 = tmp5_0*w9;
1614                              const register double tmp57_1 = tmp18_0*w2;                              const double tmp37_1 = A_10_1*w4;
1615                              const register double tmp28_1 = tmp9_0*w7;                              const double tmp38_1 = tmp5_0*w15;
1616                              const register double tmp29_1 = tmp1_0*w14;                              const double tmp17_1 = A_01_1*w7;
1617                              const register double tmp51_1 = tmp11_0*w16;                              const double tmp12_1 = tmp7_0*w4;
1618                              const register double tmp42_1 = tmp12_0*w16;                              const double tmp22_1 = tmp10_0*w7;
1619                              const register double tmp49_1 = tmp12_0*w17;                              const double tmp57_1 = tmp18_0*w2;
1620                              const register double tmp21_1 = tmp1_0*w11;                              const double tmp28_1 = tmp9_0*w7;
1621                              const register double tmp1_1 = tmp0_0*w1;                              const double tmp29_1 = tmp1_0*w14;
1622                              const register double tmp45_1 = tmp6_0*w4;                              const double tmp51_1 = tmp11_0*w16;
1623                              const register double tmp7_1 = A_10_0*w2;                              const double tmp42_1 = tmp12_0*w16;
1624                              const register double tmp6_1 = tmp3_0*w6;                              const double tmp49_1 = tmp12_0*w17;
1625                              const register double tmp13_1 = tmp8_0*w1;                              const double tmp21_1 = tmp1_0*w11;
1626                              const register double tmp36_1 = tmp16_0*w1;                              const double tmp45_1 = tmp6_0*w4;
1627                              const register double tmp41_1 = A_01_3*w2;                              const double tmp13_1 = tmp8_0*w1;
1628                              const register double tmp30_1 = tmp12_0*w13;                              const double tmp36_1 = tmp16_0*w1;
1629                              const register double tmp4_1 = A_01_2*w7;                              const double tmp41_1 = A_01_3*w2;
1630                              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;  
1631                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1632                              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;
1633                                EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
1634                              EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;                              EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
1635                              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;
1636                                EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
1637                              EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;                              EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
1638                                EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
1639                              EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1640                              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;  
1641                              EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;                              EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
1642                              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;
1643                              EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;                              EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
1644                              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;
1645                          } 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;
1646                              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;
1647                              const register double A_01 = A_p[INDEX2(0,1,2)];                          } else { // constant data
1648                              const register double A_10 = A_p[INDEX2(1,0,2)];                              const double A_00 = A_p[INDEX2(0,0,2)];
1649                              const register double A_11 = A_p[INDEX2(1,1,2)];                              const double A_10 = A_p[INDEX2(1,0,2)];
1650                              const register double tmp0_0 = A_01 + A_10;                              const double A_01 = A_p[INDEX2(0,1,2)];
1651                              const register double tmp0_1 = A_00*w18;                              const double A_11 = A_p[INDEX2(1,1,2)];
1652                              const register double tmp10_1 = A_01*w20;                              const double tmp0_0 = A_01 + A_10;
1653                              const register double tmp12_1 = A_00*w26;                              const double tmp0_1 = A_00*w18;
1654                              const register double tmp4_1 = A_00*w22;                              const double tmp1_1 = A_01*w19;
1655                              const register double tmp8_1 = A_00*w24;                              const double tmp2_1 = A_10*w20;
1656                              const register double tmp13_1 = A_10*w19;                              const double tmp3_1 = A_11*w21;
1657                              const register double tmp9_1 = tmp0_0*w20;                              const double tmp4_1 = A_00*w22;
1658                              const register double tmp3_1 = A_11*w21;                              const double tmp5_1 = tmp0_0*w19;
1659                              const register double tmp11_1 = A_11*w27;                              const double tmp6_1 = A_11*w23;
1660                              const register double tmp1_1 = A_01*w19;                              const double tmp7_1 = A_11*w25;
1661                              const register double tmp6_1 = A_11*w23;                              const double tmp8_1 = A_00*w24;
1662                              const register double tmp7_1 = A_11*w25;                              const double tmp9_1 = tmp0_0*w20;
1663                              const register double tmp2_1 = A_10*w20;                              const double tmp10_1 = A_01*w20;
1664                              const register double tmp5_1 = tmp0_0*w19;                              const double tmp11_1 = A_11*w27;
1665                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              const double tmp12_1 = A_00*w26;
1666                              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;  
1667                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1668                              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;
1669                                EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1670                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1671                              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;
1672                                EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1673                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1674                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1675                              EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;                              EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1676                              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;  
1677                              EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1678                              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;
1679                              EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1680                              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;
1681                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1682                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1683                          }                          }
1684                      }                      }
1685                      ///////////////                      ///////////////
# Line 1701  void Rectangle::assemblePDESingle(Paso_S Line 1689  void Rectangle::assemblePDESingle(Paso_S
1689                          add_EM_S=true;                          add_EM_S=true;
1690                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1691                          if (B.actsExpanded()) {                          if (B.actsExpanded()) {
1692                              const register double B_0_0 = B_p[INDEX2(0,0,2)];                              const double B_0_0 = B_p[INDEX2(0,0,2)];
1693                              const register double B_1_0 = B_p[INDEX2(1,0,2)];                              const double B_1_0 = B_p[INDEX2(1,0,2)];
1694                              const register double B_0_1 = B_p[INDEX2(0,1,2)];                              const double B_0_1 = B_p[INDEX2(0,1,2)];
1695                              const register double B_1_1 = B_p[INDEX2(1,1,2)];                              const double B_1_1 = B_p[INDEX2(1,1,2)];
1696                              const register double B_0_2 = B_p[INDEX2(0,2,2)];                              const double B_0_2 = B_p[INDEX2(0,2,2)];
1697                              const register double B_1_2 = B_p[INDEX2(1,2,2)];                              const double B_1_2 = B_p[INDEX2(1,2,2)];
1698                              const register double B_0_3 = B_p[INDEX2(0,3,2)];                              const double B_0_3 = B_p[INDEX2(0,3,2)];
1699                              const register double B_1_3 = B_p[INDEX2(1,3,2)];                              const double B_1_3 = B_p[INDEX2(1,3,2)];
1700                              const register double tmp3_0 = B_0_0 + B_0_2;                              const double tmp0_0 = B_1_0 + B_1_1;
1701                              const register double tmp1_0 = B_1_2 + B_1_3;                              const double tmp1_0 = B_1_2 + B_1_3;
1702                              const register double tmp2_0 = B_0_1 + B_0_3;                              const double tmp2_0 = B_0_1 + B_0_3;
1703                              const register double tmp0_0 = B_1_0 + B_1_1;                              const double tmp3_0 = B_0_0 + B_0_2;
1704                              const register double tmp63_1 = B_1_1*w42;                              const double tmp63_1 = B_1_1*w42;
1705                              const register double tmp79_1 = B_1_1*w40;                              const double tmp79_1 = B_1_1*w40;
1706                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
1707                              const register double tmp8_1 = tmp0_0*w32;                              const double tmp8_1 = tmp0_0*w32;
1708                              const register double tmp71_1 = B_0_1*w34;                              const double tmp71_1 = B_0_1*w34;
1709                              const register double tmp19_1 = B_0_3*w31;                              const double tmp19_1 = B_0_3*w31;
1710                              const register double tmp15_1 = B_0_3*w34;                              const double tmp15_1 = B_0_3*w34;
1711                              const register double tmp9_1 = tmp3_0*w34;                              const double tmp9_1 = tmp3_0*w34;
1712                              const register double tmp35_1 = B_1_0*w36;                              const double tmp35_1 = B_1_0*w36;
1713                              const register double tmp66_1 = B_0_3*w28;                              const double tmp66_1 = B_0_3*w28;
1714                              const register double tmp28_1 = B_1_0*w42;                              const double tmp28_1 = B_1_0*w42;
1715                              const register double tmp22_1 = B_1_0*w40;                              const double tmp22_1 = B_1_0*w40;
1716                              const register double tmp16_1 = B_1_2*w29;                              const double tmp16_1 = B_1_2*w29;
1717                              const register double tmp6_1 = tmp2_0*w35;                              const double tmp6_1 = tmp2_0*w35;
1718                              const register double tmp55_1 = B_1_3*w40;                              const double tmp55_1 = B_1_3*w40;
1719                              const register double tmp50_1 = B_1_3*w42;                              const double tmp50_1 = B_1_3*w42;
1720                              const register double tmp7_1 = tmp1_0*w29;                              const double tmp7_1 = tmp1_0*w29;
1721                              const register double tmp1_1 = tmp1_0*w32;                              const double tmp1_1 = tmp1_0*w32;
1722                              const register double tmp57_1 = B_0_3*w30;                              const double tmp57_1 = B_0_3*w30;
1723                              const register double tmp18_1 = B_1_1*w32;                              const double tmp18_1 = B_1_1*w32;
1724                              const register double tmp53_1 = B_1_0*w41;                              const double tmp53_1 = B_1_0*w41;
1725                              const register double tmp61_1 = B_1_3*w36;                              const double tmp61_1 = B_1_3*w36;
1726                              const register double tmp27_1 = B_0_3*w38;                              const double tmp27_1 = B_0_3*w38;
1727                              const register double tmp64_1 = B_0_2*w30;                              const double tmp64_1 = B_0_2*w30;
1728                              const register double tmp76_1 = B_0_1*w38;                              const double tmp76_1 = B_0_1*w38;
1729                              const register double tmp39_1 = tmp2_0*w34;                              const double tmp39_1 = tmp2_0*w34;
1730                              const register double tmp62_1 = B_0_1*w31;                              const double tmp62_1 = B_0_1*w31;
1731                              const register double tmp56_1 = B_0_0*w31;                              const double tmp56_1 = B_0_0*w31;
1732                              const register double tmp49_1 = B_1_1*w36;                              const double tmp49_1 = B_1_1*w36;
1733                              const register double tmp2_1 = B_0_2*w31;                              const double tmp2_1 = B_0_2*w31;
1734                              const register double tmp23_1 = B_0_2*w33;                              const double tmp23_1 = B_0_2*w33;
1735                              const register double tmp38_1 = B_1_1*w43;                              const double tmp38_1 = B_1_1*w43;
1736                              const register double tmp74_1 = B_1_2*w41;                              const double tmp74_1 = B_1_2*w41;
1737                              const register double tmp43_1 = B_1_1*w41;                              const double tmp43_1 = B_1_1*w41;
1738                              const register double tmp58_1 = B_0_2*w28;                              const double tmp58_1 = B_0_2*w28;
1739                              const register double tmp67_1 = B_0_0*w33;                              const double tmp67_1 = B_0_0*w33;
1740                              const register double tmp33_1 = tmp0_0*w39;                              const double tmp33_1 = tmp0_0*w39;
1741                              const register double tmp4_1 = B_0_0*w28;                              const double tmp4_1 = B_0_0*w28;
1742                              const register double tmp20_1 = B_0_0*w30;                              const double tmp20_1 = B_0_0*w30;
1743                              const register double tmp13_1 = B_0_2*w38;                              const double tmp13_1 = B_0_2*w38;
1744                              const register double tmp65_1 = B_1_2*w43;                              const double tmp65_1 = B_1_2*w43;
1745                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
1746                              const register double tmp41_1 = tmp3_0*w33;                              const double tmp41_1 = tmp3_0*w33;
1747                              const register double tmp73_1 = B_0_2*w37;                              const double tmp73_1 = B_0_2*w37;
1748                              const register double tmp69_1 = B_0_0*w38;                              const double tmp69_1 = B_0_0*w38;
1749                              const register double tmp48_1 = B_1_2*w39;                              const double tmp48_1 = B_1_2*w39;
1750                              const register double tmp59_1 = B_0_1*w33;                              const double tmp59_1 = B_0_1*w33;
1751                              const register double tmp17_1 = B_1_3*w41;                              const double tmp17_1 = B_1_3*w41;
1752                              const register double tmp5_1 = B_0_3*w33;                              const double tmp5_1 = B_0_3*w33;
1753                              const register double tmp3_1 = B_0_1*w30;                              const double tmp3_1 = B_0_1*w30;
1754                              const register double tmp21_1 = B_0_1*w28;                              const double tmp21_1 = B_0_1*w28;
1755                              const register double tmp42_1 = B_1_0*w29;                              const double tmp42_1 = B_1_0*w29;
1756                              const register double tmp54_1 = B_1_2*w32;                              const double tmp54_1 = B_1_2*w32;
1757                              const register double tmp60_1 = B_1_0*w39;                              const double tmp60_1 = B_1_0*w39;
1758                              const register double tmp32_1 = tmp1_0*w36;                              const double tmp32_1 = tmp1_0*w36;
1759                              const register double tmp10_1 = B_0_1*w37;                              const double tmp10_1 = B_0_1*w37;
1760                              const register double tmp14_1 = B_0_0*w35;                              const double tmp14_1 = B_0_0*w35;
1761                              const register double tmp29_1 = B_0_1*w35;                              const double tmp29_1 = B_0_1*w35;
1762                              const register double tmp26_1 = B_1_2*w36;                              const double tmp26_1 = B_1_2*w36;
1763                              const register double tmp30_1 = B_1_3*w43;                              const double tmp30_1 = B_1_3*w43;
1764                              const register double tmp70_1 = B_0_2*w35;                              const double tmp70_1 = B_0_2*w35;
1765                              const register double tmp34_1 = B_1_3*w39;                              const double tmp34_1 = B_1_3*w39;
1766                              const register double tmp51_1 = B_1_0*w43;                              const double tmp51_1 = B_1_0*w43;
1767                              const register double tmp31_1 = B_0_2*w34;                              const double tmp31_1 = B_0_2*w34;
1768                              const register double tmp45_1 = tmp3_0*w28;                              const double tmp45_1 = tmp3_0*w28;
1769                              const register double tmp11_1 = tmp1_0*w39;                              const double tmp11_1 = tmp1_0*w39;
1770                              const register double tmp52_1 = B_1_1*w29;                              const double tmp52_1 = B_1_1*w29;
1771                              const register double tmp44_1 = B_1_3*w32;                              const double tmp44_1 = B_1_3*w32;
1772                              const register double tmp25_1 = B_1_1*w39;                              const double tmp25_1 = B_1_1*w39;
1773                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
1774                              const register double tmp72_1 = B_1_3*w29;                              const double tmp72_1 = B_1_3*w29;
1775                              const register double tmp40_1 = tmp2_0*w28;                              const double tmp40_1 = tmp2_0*w28;
1776                              const register double tmp46_1 = B_1_2*w40;                              const double tmp46_1 = B_1_2*w40;
1777                              const register double tmp36_1 = B_1_2*w42;                              const double tmp36_1 = B_1_2*w42;
1778                              const register double tmp24_1 = B_0_0*w37;                              const double tmp24_1 = B_0_0*w37;
1779                              const register double tmp77_1 = B_0_3*w35;                              const double tmp77_1 = B_0_3*w35;
1780                              const register double tmp68_1 = B_0_3*w37;                              const double tmp68_1 = B_0_3*w37;
1781                              const register double tmp78_1 = B_0_0*w34;                              const double tmp78_1 = B_0_0*w34;
1782                              const register double tmp12_1 = tmp0_0*w36;                              const double tmp12_1 = tmp0_0*w36;
1783                              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;  
1784                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1785                              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;
1786                                EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1787                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
1788                              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;
1789                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1790                              EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;                              EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1791                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1792                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1793                              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;  
1794                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1795                              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;
1796                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
1797                              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;
1798                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1799                              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;
1800                              const register double B_1 = B_p[1];                          } else { // constant data
1801                              const register double tmp6_1 = B_1*w50;                              const double B_0 = B_p[0];
1802                              const register double tmp1_1 = B_1*w45;                              const double B_1 = B_p[1];
1803                              const register double tmp5_1 = B_1*w49;                              const double tmp0_1 = B_0*w44;
1804                              const register double tmp4_1 = B_1*w48;                              const double tmp1_1 = B_1*w45;
1805                              const register double tmp0_1 = B_0*w44;                              const double tmp2_1 = B_0*w46;
1806                              const register double tmp2_1 = B_0*w46;                              const double tmp3_1 = B_0*w47;
1807                              const register double tmp7_1 = B_0*w51;                              const double tmp4_1 = B_1*w48;
1808                              const register double tmp3_1 = B_0*w47;                              const double tmp5_1 = B_1*w49;
1809                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const double tmp6_1 = B_1*w50;
1810                              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;  
1811                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1812                              EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1813                                EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
1814                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
1815                              EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1816                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
1817                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
1818                                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
1819                              EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;                              EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
1820                              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;  
1821                              EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;                              EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
1822                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1823                              EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;                              EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
1824                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
1825                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
1826                                EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
1827                          }                          }
1828                      }                      }
1829                      ///////////////                      ///////////////
# Line 1845  void Rectangle::assemblePDESingle(Paso_S Line 1833  void Rectangle::assemblePDESingle(Paso_S
1833                          add_EM_S=true;                          add_EM_S=true;
1834                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
1835                          if (C.actsExpanded()) {                          if (C.actsExpanded()) {
1836                              const register double C_0_0 = C_p[INDEX2(0,0,2)];                              const double C_0_0 = C_p[INDEX2(0,0,2)];
1837                              const register double C_1_0 = C_p[INDEX2(1,0,2)];                              const double C_1_0 = C_p[INDEX2(1,0,2)];
1838                              const register double C_0_1 = C_p[INDEX2(0,1,2)];                              const double C_0_1 = C_p[INDEX2(0,1,2)];
1839                              const register double C_1_1 = C_p[INDEX2(1,1,2)];                              const double C_1_1 = C_p[INDEX2(1,1,2)];
1840                              const register double C_0_2 = C_p[INDEX2(0,2,2)];                              const double C_0_2 = C_p[INDEX2(0,2,2)];
1841                              const register double C_1_2 = C_p[INDEX2(1,2,2)];                              const double C_1_2 = C_p[INDEX2(1,2,2)];
1842                              const register double C_0_3 = C_p[INDEX2(0,3,2)];                              const double C_0_3 = C_p[INDEX2(0,3,2)];
1843                              const register double C_1_3 = C_p[INDEX2(1,3,2)];                              const double C_1_3 = C_p[INDEX2(1,3,2)];
1844                              const register double tmp2_0 = C_0_1 + C_0_3;                              const double tmp0_0 = C_1_0 + C_1_1;
1845                              const register double tmp1_0 = C_1_2 + C_1_3;                              const double tmp1_0 = C_1_2 + C_1_3;
1846                              const register double tmp3_0 = C_0_0 + C_0_2;                              const double tmp2_0 = C_0_1 + C_0_3;
1847                              const register double tmp0_0 = C_1_0 + C_1_1;                              const double tmp3_0 = C_0_0 + C_0_2;
1848                              const register double tmp64_1 = C_0_2*w30;                              const double tmp64_1 = C_0_2*w30;
1849                              const register double tmp14_1 = C_0_2*w28;                              const double tmp14_1 = C_0_2*w28;
1850                              const register double tmp19_1 = C_0_3*w31;                              const double tmp19_1 = C_0_3*w31;
1851                              const register double tmp22_1 = C_1_0*w40;                              const double tmp22_1 = C_1_0*w40;
1852                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
1853                              const register double tmp29_1 = C_0_1*w35;                              const double tmp29_1 = C_0_1*w35;
1854                              const register double tmp73_1 = C_0_2*w37;                              const double tmp73_1 = C_0_2*w37;
1855                              const register double tmp74_1 = C_1_2*w41;                              const double tmp74_1 = C_1_2*w41;
1856                              const register double tmp52_1 = C_1_3*w39;                              const double tmp52_1 = C_1_3*w39;
1857                              const register double tmp25_1 = C_1_1*w39;                              const double tmp25_1 = C_1_1*w39;
1858                              const register double tmp62_1 = C_0_1*w31;                              const double tmp62_1 = C_0_1*w31;
1859                              const register double tmp79_1 = C_1_1*w40;                              const double tmp79_1 = C_1_1*w40;
1860                              const register double tmp43_1 = C_1_1*w36;                              const double tmp43_1 = C_1_1*w36;
1861                              const register double tmp27_1 = C_0_3*w38;                              const double tmp27_1 = C_0_3*w38;
1862                              const register double tmp28_1 = C_1_0*w42;                              const double tmp28_1 = C_1_0*w42;
1863                              const register double tmp63_1 = C_1_1*w42;                              const double tmp63_1 = C_1_1*w42;
1864                              const register double tmp59_1 = C_0_3*w34;                              const double tmp59_1 = C_0_3*w34;
1865                              const register double tmp72_1 = C_1_3*w29;                              const double tmp72_1 = C_1_3*w29;
1866                              const register double tmp40_1 = tmp2_0*w35;                              const double tmp40_1 = tmp2_0*w35;
1867                              const register double tmp13_1 = C_0_3*w30;                              const double tmp13_1 = C_0_3*w30;
1868                              const register double tmp51_1 = C_1_2*w40;                              const double tmp51_1 = C_1_2*w40;
1869                              const register double tmp54_1 = C_1_2*w42;                              const double tmp54_1 = C_1_2*w42;
1870                              const register double tmp12_1 = C_0_0*w31;                              const double tmp12_1 = C_0_0*w31;
1871                              const register double tmp2_1 = tmp1_0*w32;                              const double tmp2_1 = tmp1_0*w32;
1872                              const register double tmp68_1 = C_0_2*w31;                              const double tmp68_1 = C_0_2*w31;
1873                              const register double tmp75_1 = C_1_0*w32;                              const double tmp75_1 = C_1_0*w32;
1874                              const register double tmp49_1 = C_1_1*w41;                              const double tmp49_1 = C_1_1*w41;
1875                              const register double tmp4_1 = C_0_2*w35;                              const double tmp4_1 = C_0_2*w35;
1876                              const register double tmp66_1 = C_0_3*w28;                              const double tmp66_1 = C_0_3*w28;
1877                              const register double tmp56_1 = C_0_1*w37;                              const double tmp56_1 = C_0_1*w37;
1878                              const register double tmp5_1 = C_0_1*w34;                              const double tmp5_1 = C_0_1*w34;
1879                              const register double tmp38_1 = tmp2_0*w34;                              const double tmp38_1 = tmp2_0*w34;
1880                              const register double tmp76_1 = C_0_1*w38;                              const double tmp76_1 = C_0_1*w38;
1881                              const register double tmp21_1 = C_0_1*w28;                              const double tmp21_1 = C_0_1*w28;
1882                              const register double tmp69_1 = C_0_1*w30;                              const double tmp69_1 = C_0_1*w30;
1883                              const register double tmp53_1 = C_1_0*w36;                              const double tmp53_1 = C_1_0*w36;
1884                              const register double tmp42_1 = C_1_2*w39;                              const double tmp42_1 = C_1_2*w39;
1885                              const register double tmp32_1 = tmp1_0*w29;                              const double tmp32_1 = tmp1_0*w29;
1886                              const register double tmp45_1 = C_1_0*w43;                              const double tmp45_1 = C_1_0*w43;
1887                              const register double tmp33_1 = tmp0_0*w32;                              const double tmp33_1 = tmp0_0*w32;
1888                              const register double tmp35_1 = C_1_0*w41;                              const double tmp35_1 = C_1_0*w41;
1889                              const register double tmp26_1 = C_1_2*w36;                              const double tmp26_1 = C_1_2*w36;
1890                              const register double tmp67_1 = C_0_0*w33;                              const double tmp67_1 = C_0_0*w33;
1891                              const register double tmp31_1 = C_0_2*w34;                              const double tmp31_1 = C_0_2*w34;
1892                              const register double tmp20_1 = C_0_0*w30;                              const double tmp20_1 = C_0_0*w30;
1893                              const register double tmp70_1 = C_0_0*w28;                              const double tmp70_1 = C_0_0*w28;
1894                              const register double tmp8_1 = tmp0_0*w39;                              const double tmp8_1 = tmp0_0*w39;
1895                              const register double tmp30_1 = C_1_3*w43;                              const double tmp30_1 = C_1_3*w43;
1896                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
1897                              const register double tmp17_1 = C_1_3*w41;                              const double tmp17_1 = C_1_3*w41;
1898                              const register double tmp58_1 = C_0_0*w35;                              const double tmp58_1 = C_0_0*w35;
1899                              const register double tmp9_1 = tmp3_0*w33;                              const double tmp9_1 = tmp3_0*w33;
1900                              const register double tmp61_1 = C_1_3*w36;                              const double tmp61_1 = C_1_3*w36;
1901                              const register double tmp41_1 = tmp3_0*w34;                              const double tmp41_1 = tmp3_0*w34;
1902                              const register double tmp50_1 = C_1_3*w32;                              const double tmp50_1 = C_1_3*w32;
1903                              const register double tmp18_1 = C_1_1*w32;                              const double tmp18_1 = C_1_1*w32;
1904                              const register double tmp6_1 = tmp1_0*w36;                              const double tmp6_1 = tmp1_0*w36;
1905                              const register double tmp3_1 = C_0_0*w38;                              const double tmp3_1 = C_0_0*w38;
1906                              const register double tmp34_1 = C_1_1*w29;                              const double tmp34_1 = C_1_1*w29;
1907                              const register double tmp77_1 = C_0_3*w35;                              const double tmp77_1 = C_0_3*w35;
1908                              const register double tmp65_1 = C_1_2*w43;                              const double tmp65_1 = C_1_2*w43;
1909                              const register double tmp71_1 = C_0_3*w33;                              const double tmp71_1 = C_0_3*w33;
1910                              const register double tmp55_1 = C_1_1*w43;                              const double tmp55_1 = C_1_1*w43;
1911                              const register double tmp46_1 = tmp3_0*w28;                              const double tmp46_1 = tmp3_0*w28;
1912                              const register double tmp24_1 = C_0_0*w37;                              const double tmp24_1 = C_0_0*w37;
1913                              const register double tmp10_1 = tmp1_0*w39;                              const double tmp10_1 = tmp1_0*w39;
1914                              const register double tmp48_1 = C_1_0*w29;                              const double tmp48_1 = C_1_0*w29;
1915                              const register double tmp15_1 = C_0_1*w33;                              const double tmp15_1 = C_0_1*w33;
1916                              const register double tmp36_1 = C_1_2*w32;                              const double tmp36_1 = C_1_2*w32;
1917                              const register double tmp60_1 = C_1_0*w39;                              const double tmp60_1 = C_1_0*w39;
1918                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
1919                              const register double tmp16_1 = C_1_2*w29;                              const double tmp16_1 = C_1_2*w29;
1920                              const register double tmp1_1 = C_0_3*w37;                              const double tmp1_1 = C_0_3*w37;
1921                              const register double tmp7_1 = tmp2_0*w28;                              const double tmp7_1 = tmp2_0*w28;
1922                              const register double tmp39_1 = C_1_3*w40;                              const double tmp39_1 = C_1_3*w40;
1923                              const register double tmp44_1 = C_1_3*w42;                              const double tmp44_1 = C_1_3*w42;
1924                              const register double tmp57_1 = C_0_2*w38;                              const double tmp57_1 = C_0_2*w38;
1925                              const register double tmp78_1 = C_0_0*w34;                              const double tmp78_1 = C_0_0*w34;
1926                              const register double tmp11_1 = tmp0_0*w36;                              const double tmp11_1 = tmp0_0*w36;
1927                              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;  
1928                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1929                              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;
1930                                EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1931                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
1932                              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;
1933                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1934                              EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;                              EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1935                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1936                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1937                              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;  
1938                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1939                              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;
1940                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
1941                              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;
1942                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1943                              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;
1944                              const register double C_1 = C_p[1];                          } else { // constant data
1945                              const register double tmp1_1 = C_1*w45;                              const double C_0 = C_p[0];
1946                              const register double tmp3_1 = C_0*w51;                              const double C_1 = C_p[1];
1947                              const register double tmp4_1 = C_0*w44;                              const double tmp0_1 = C_0*w47;
1948                              const register double tmp7_1 = C_0*w46;                              const double tmp1_1 = C_1*w45;
1949                              const register double tmp5_1 = C_1*w49;                              const double tmp2_1 = C_1*w48;
1950                              const register double tmp2_1 = C_1*w48;                              const double tmp3_1 = C_0*w51;
1951                              const register double tmp0_1 = C_0*w47;                              const double tmp4_1 = C_0*w44;
1952                              const register double tmp6_1 = C_1*w50;                              const double tmp5_1 = C_1*w49;
1953                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const double tmp6_1 = C_1*w50;
1954                              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;  
1955                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1956                              EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
1957                                EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
1958                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
1959                              EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1960                                EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
1961                              EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;                              EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
1962                                EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
1963                              EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;                              EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
1964                              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;  
1965                              EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;                              EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
1966                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;                              EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
1967                              EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;                              EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
1968                              EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
1969                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
1970                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
1971                          }                          }
1972                      }                      }
1973                      ///////////////                      ///////////////
# Line 1989  void Rectangle::assemblePDESingle(Paso_S Line 1977  void Rectangle::assemblePDESingle(Paso_S
1977                          add_EM_S=true;                          add_EM_S=true;
1978                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
1979                          if (D.actsExpanded()) {                          if (D.actsExpanded()) {
1980                              const register double D_0 = D_p[0];                              const double D_0 = D_p[0];
1981                              const register double D_1 = D_p[1];                              const double D_1 = D_p[1];
1982                              const register double D_2 = D_p[2];                              const double D_2 = D_p[2];
1983                              const register double D_3 = D_p[3];                              const double D_3 = D_p[3];
1984                              const register double tmp4_0 = D_1 + D_3;                              const double tmp0_0 = D_0 + D_1;
1985                              const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;                              const double tmp1_0 = D_2 + D_3;
1986                              const register double tmp5_0 = D_0 + D_2;                              const double tmp2_0 = D_0 + D_1 + D_2 + D_3;
1987                              const register double tmp0_0 = D_0 + D_1;                              const double tmp3_0 = D_1 + D_2;
1988                              const register double tmp6_0 = D_0 + D_3;                              const double tmp4_0 = D_1 + D_3;
1989                              const register double tmp1_0 = D_2 + D_3;                              const double tmp5_0 = D_0 + D_2;
1990                              const register double tmp3_0 = D_1 + D_2;                              const double tmp6_0 = D_0 + D_3;
1991                              const register double tmp16_1 = D_1*w56;                              const double tmp0_1 = tmp0_0*w52;
1992                              const register double tmp14_1 = tmp6_0*w54;                              const double tmp1_1 = tmp1_0*w53;
1993                              const register double tmp8_1 = D_3*w55;                              const double tmp2_1 = tmp2_0*w54;
1994                              const register double tmp2_1 = tmp2_0*w54;                              const double tmp3_1 = tmp1_0*w52;
1995                              const register double tmp12_1 = tmp5_0*w52;                              const double tmp4_1 = tmp0_0*w53;
1996                              const register double tmp4_1 = tmp0_0*w53;                              const double tmp5_1 = tmp3_0*w54;
1997                              const register double tmp3_1 = tmp1_0*w52;                              const double tmp6_1 = D_0*w55;
1998                              const register double tmp13_1 = tmp4_0*w53;                              const double tmp7_1 = D_3*w56;
1999                              const register double tmp10_1 = tmp4_0*w52;                              const double tmp8_1 = D_3*w55;
2000                              const register double tmp1_1 = tmp1_0*w53;                              const double tmp9_1 = D_0*w56;
2001                              const register double tmp7_1 = D_3*w56;                              const double tmp10_1 = tmp4_0*w52;
2002                              const register double tmp0_1 = tmp0_0*w52;                              const double tmp11_1 = tmp5_0*w53;
2003                              const register double tmp11_1 = tmp5_0*w53;                              const double tmp12_1 = tmp5_0*w52;
2004                              const register double tmp9_1 = D_0*w56;                              const double tmp13_1 = tmp4_0*w53;
2005                              const register double tmp5_1 = tmp3_0*w54;                              const double tmp14_1 = tmp6_0*w54;
2006                              const register double tmp18_1 = D_2*w56;                              const double tmp15_1 = D_2*w55;
2007                              const register double tmp17_1 = D_1*w55;                              const double tmp16_1 = D_1*w56;
2008                              const register double tmp6_1 = D_0*w55;                              const double tmp17_1 = D_1*w55;
2009                              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;  
2010                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2011                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2012                                EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
2013                              EM_S[INDEX2(3,0,4)]+=tmp2_1;                              EM_S[INDEX2(3,0,4)]+=tmp2_1;
2014                              EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2015                                EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
2016                              EM_S[INDEX2(2,1,4)]+=tmp2_1;                              EM_S[INDEX2(2,1,4)]+=tmp2_1;
2017                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
2018                              EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;                              EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
2019                              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;  
2020                              EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;                              EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
2021                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
2022                              EM_S[INDEX2(0,3,4)]+=tmp2_1;                              EM_S[INDEX2(0,3,4)]+=tmp2_1;
2023                              EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;                              EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
2024                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2025                              const register double D_0 = D_p[0];                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2026                              const register double tmp2_1 = D_0*w59;                          } else { // constant data
2027                              const register double tmp1_1 = D_0*w58;                              const double tmp0_1 = D_p[0]*w57;
2028                              const register double tmp0_1 = D_0*w57;                              const double tmp1_1 = D_p[0]*w58;
2029                              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;  
2030                              EM_S[INDEX2(0,0,4)]+=tmp2_1;                              EM_S[INDEX2(0,0,4)]+=tmp2_1;
2031                              EM_S[INDEX2(3,3,4)]+=tmp2_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
2032                                EM_S[INDEX2(2,0,4)]+=tmp0_1;
2033                              EM_S[INDEX2(3,0,4)]+=tmp1_1;                              EM_S[INDEX2(3,0,4)]+=tmp1_1;
2034                              EM_S[INDEX2(3,1,4)]+=tmp0_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1;
2035                                EM_S[INDEX2(1,1,4)]+=tmp2_1;
2036                              EM_S[INDEX2(2,1,4)]+=tmp1_1;                              EM_S[INDEX2(2,1,4)]+=tmp1_1;
2037                                EM_S[INDEX2(3,1,4)]+=tmp0_1;
2038                              EM_S[INDEX2(0,2,4)]+=tmp0_1;                              EM_S[INDEX2(0,2,4)]+=tmp0_1;
2039                              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;  
2040                              EM_S[INDEX2(2,2,4)]+=tmp2_1;                              EM_S[INDEX2(2,2,4)]+=tmp2_1;
2041                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(3,2,4)]+=tmp0_1;
2042                              EM_S[INDEX2(0,3,4)]+=tmp1_1;                              EM_S[INDEX2(0,3,4)]+=tmp1_1;
2043                              EM_S[INDEX2(1,1,4)]+=tmp2_1;                              EM_S[INDEX2(1,3,4)]+=tmp0_1;
2044                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
2045                                EM_S[INDEX2(3,3,4)]+=tmp2_1;
2046                          }                          }
2047                      }                      }
2048                      ///////////////                      ///////////////
# Line 2065  void Rectangle::assemblePDESingle(Paso_S Line 2052  void Rectangle::assemblePDESingle(Paso_S
2052                          add_EM_F=true;                          add_EM_F=true;
2053                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2054                          if (X.actsExpanded()) {                          if (X.actsExpanded()) {
2055                              const register double X_0_0 = X_p[INDEX2(0,0,2)];                              const double X_0_0 = X_p[INDEX2(0,0,2)];
2056                              const register double X_1_0 = X_p[INDEX2(1,0,2)];                              const double X_1_0 = X_p[INDEX2(1,0,2)];
2057                              const register double X_0_1 = X_p[INDEX2(0,1,2)];                              const double X_0_1 = X_p[INDEX2(0,1,2)];
2058                              const register double X_1_1 = X_p[INDEX2(1,1,2)];                              const double X_1_1 = X_p[INDEX2(1,1,2)];
2059                              const register double X_0_2 = X_p[INDEX2(0,2,2)];                              const double X_0_2 = X_p[INDEX2(0,2,2)];
2060                              const register double X_1_2 = X_p[INDEX2(1,2,2)];                              const double X_1_2 = X_p[INDEX2(1,2,2)];
2061                              const register double X_0_3 = X_p[INDEX2(0,3,2)];                              const double X_0_3 = X_p[INDEX2(0,3,2)];
2062                              const register double X_1_3 = X_p[INDEX2(1,3,2)];                              const double X_1_3 = X_p[INDEX2(1,3,2)];
2063                              const register double tmp1_0 = X_1_1 + X_1_3;                              const double tmp0_0 = X_0_2 + X_0_3;
2064                              const register double tmp3_0 = X_0_0 + X_0_1;                              const double tmp1_0 = X_1_1 + X_1_3;
2065                              const register double tmp2_0 = X_1_0 + X_1_2;                              const double tmp2_0 = X_1_0 + X_1_2;
2066                              const register double tmp0_0 = X_0_2 + X_0_3;                              const double tmp3_0 = X_0_0 + X_0_1;
2067                              const register double tmp8_1 = tmp2_0*w66;                              const double tmp0_1 = tmp0_0*w63;
2068                              const register double tmp5_1 = tmp3_0*w64;                              const double tmp1_1 = tmp1_0*w62;
2069                              const register double tmp14_1 = tmp0_0*w64;                              const double tmp2_1 = tmp2_0*w61;
2070                              const register double tmp3_1 = tmp3_0*w60;                              const double tmp3_1 = tmp3_0*w60;
2071                              const register double tmp9_1 = tmp3_0*w63;                              const double tmp4_1 = tmp0_0*w65;
2072                              const register double tmp13_1 = tmp3_0*w65;                              const double tmp5_1 = tmp3_0*w64;
2073                              const register double tmp12_1 = tmp1_0*w66;                              const double tmp6_1 = tmp2_0*w62;
2074                              const register double tmp10_1 = tmp0_0*w60;                              const double tmp7_1 = tmp1_0*w61;
2075                              const register double tmp2_1 = tmp2_0*w61;                              const double tmp8_1 = tmp2_0*w66;
2076                              const register double tmp6_1 = tmp2_0*w62;                              const double tmp9_1 = tmp3_0*w63;
2077                              const register double tmp4_1 = tmp0_0*w65;                              const double tmp10_1 = tmp0_0*w60;
2078                              const register double tmp11_1 = tmp1_0*w67;                              const double tmp11_1 = tmp1_0*w67;
2079                              const register double tmp1_1 = tmp1_0*w62;                              const double tmp12_1 = tmp1_0*w66;
2080                              const register double tmp7_1 = tmp1_0*w61;                              const double tmp13_1 = tmp3_0*w65;
2081                              const register double tmp0_1 = tmp0_0*w63;                              const double tmp14_1 = tmp0_0*w64;
2082                              const register double tmp15_1 = tmp2_0*w67;                              const double tmp15_1 = tmp2_0*w67;
2083                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2084                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2085                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2086                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2087                          } else { /* constant data */                          } else { // constant data
2088                              const register double X_0 = X_p[0];                              const double tmp0_1 = X_p[1]*w69;
2089                              const register double X_1 = X_p[1];                              const double tmp1_1 = X_p[0]*w68;
2090                              const register double tmp3_1 = X_1*w71;                              const double tmp2_1 = X_p[0]*w70;
2091                              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;  
2092                              EM_F[0]+=tmp0_1 + tmp1_1;                              EM_F[0]+=tmp0_1 + tmp1_1;
2093                              EM_F[1]+=tmp0_1 + tmp2_1;                              EM_F[1]+=tmp0_1 + tmp2_1;
2094                              EM_F[2]+=tmp1_1 + tmp3_1;                              EM_F[2]+=tmp1_1 + tmp3_1;
# Line 2117  void Rectangle::assemblePDESingle(Paso_S Line 2102  void Rectangle::assemblePDESingle(Paso_S
2102                          add_EM_F=true;                          add_EM_F=true;
2103                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2104                          if (Y.actsExpanded()) {                          if (Y.actsExpanded()) {
2105                              const register double Y_0 = Y_p[0];                              const double Y_0 = Y_p[0];
2106                              const register double Y_1 = Y_p[1];                              const double Y_1 = Y_p[1];
2107                              const register double Y_2 = Y_p[2];                              const double Y_2 = Y_p[2];
2108                              const register double Y_3 = Y_p[3];                              const double Y_3 = Y_p[3];
2109                              const register double tmp0_0 = Y_1 + Y_2;                              const double tmp0_0 = Y_1 + Y_2;
2110                              const register double tmp1_0 = Y_0 + Y_3;                              const double tmp1_0 = Y_0 + Y_3;
2111                              const register double tmp9_1 = Y_0*w74;                              const double tmp0_1 = Y_0*w72;
2112                              const register double tmp4_1 = tmp1_0*w73;                              const double tmp1_1 = tmp0_0*w73;
2113                              const register double tmp5_1 = Y_2*w74;                              const double tmp2_1 = Y_3*w74;
2114                              const register double tmp7_1 = Y_1*w74;                              const double tmp3_1 = Y_1*w72;
2115                              const register double tmp6_1 = Y_2*w72;                              const double tmp4_1 = tmp1_0*w73;
2116                              const register double tmp2_1 = Y_3*w74;                              const double tmp5_1 = Y_2*w74;
2117                              const register double tmp8_1 = Y_3*w72;                              const double tmp6_1 = Y_2*w72;
2118                              const register double tmp3_1 = Y_1*w72;                              const double tmp7_1 = Y_1*w74;
2119                              const register double tmp0_1 = Y_0*w72;                              const double tmp8_1 = Y_3*w72;
2120                              const register double tmp1_1 = tmp0_0*w73;                              const double tmp9_1 = Y_0*w74;
2121                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2122                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2123                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2124                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2125                          } else { /* constant data */                          } else { // constant data
2126                              const register double Y_0 = Y_p[0];                              const double tmp0_1 = Y_p[0]*w75;
                             const register double tmp0_1 = Y_0*w75;  
2127                              EM_F[0]+=tmp0_1;                              EM_F[0]+=tmp0_1;
2128                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
2129                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
2130                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
2131                          }                          }
2132                      }                      }
                     /* GENERATOR SNIP_PDE_SINGLE BOTTOM */  
2133    
2134                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2135                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
# Line 2156  void Rectangle::assemblePDESingle(Paso_S Line 2139  void Rectangle::assemblePDESingle(Paso_S
2139                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2140                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2141                      if (add_EM_F) {                      if (add_EM_F) {
                         //cout << "-- AddtoRHS -- " << endl;  
2142                          double *F_p=rhs.getSampleDataRW(0);                          double *F_p=rhs.getSampleDataRW(0);
2143                          for (index_t i=0; i<rowIndex.size(); i++) {                          for (index_t i=0; i<rowIndex.size(); i++) {
2144                              if (rowIndex[i]<getNumDOF()) {                              if (rowIndex[i]<getNumDOF()) {
2145                                  F_p[rowIndex[i]]+=EM_F[i];                                  F_p[rowIndex[i]]+=EM_F[i];
                                 //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;  
2146                              }                              }
2147                          }                          }
                         //cout << "---"<<endl;  
2148                      }                      }
2149                      if (add_EM_S) {                      if (add_EM_S) {
                         //cout << "-- AddtoSystem -- " << endl;  
2150                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2151                      }                      }
2152    
2153                  } // end k0 loop                  } // end k0 loop
2154              } // end k1 loop              } // end k1 loop
2155          } // end of coloring          } // end of colouring
2156        } // end of parallel region
2157    }
2158    
2159    //protected
2160    void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2161            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2162            const escript::Data& C, const escript::Data& D,
2163            const escript::Data& X, const escript::Data& Y) const
2164    {
2165        const double h0 = m_l0/m_gNE0;
2166        const double h1 = m_l1/m_gNE1;
2167        const double w0 = -.25*h1/h0;
2168        const double w1 = .25;
2169        const double w2 = -.25;
2170        const double w3 = .25*h0/h1;
2171        const double w4 = -.25*h0/h1;
2172        const double w5 = .25*h1/h0;
2173        const double w6 = -.125*h1;
2174        const double w7 = -.125*h0;
2175        const double w8 = .125*h1;
2176        const double w9 = .125*h0;
2177        const double w10 = .0625*h0*h1;
2178        const double w11 = -.5*h1;
2179        const double w12 = -.5*h0;
2180        const double w13 = .5*h1;
2181        const double w14 = .5*h0;
2182        const double w15 = .25*h0*h1;
2183    
2184        rhs.requireWrite();
2185    #pragma omp parallel
2186        {
2187            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2188    #pragma omp for
2189                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2190                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2191                        bool add_EM_S=false;
2192                        bool add_EM_F=false;
2193                        vector<double> EM_S(4*4, 0);
2194                        vector<double> EM_F(4, 0);
2195                        const index_t e = k0 + m_NE0*k1;
2196                        ///////////////
2197                        // process A //
2198                        ///////////////
2199                        if (!A.isEmpty()) {
2200                            add_EM_S=true;
2201                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2202                            const double A_00 = A_p[INDEX2(0,0,2)];
2203                            const double A_10 = A_p[INDEX2(1,0,2)];
2204                            const double A_01 = A_p[INDEX2(0,1,2)];
2205                            const double A_11 = A_p[INDEX2(1,1,2)];
2206                            const double tmp0_0 = A_01 + A_10;
2207                            const double tmp0_1 = A_11*w3;
2208                            const double tmp1_1 = A_00*w0;
2209                            const double tmp2_1 = A_01*w1;
2210                            const double tmp3_1 = A_10*w2;
2211                            const double tmp4_1 = tmp0_0*w1;
2212                            const double tmp5_1 = A_11*w4;
2213                            const double tmp6_1 = A_00*w5;
2214                            const double tmp7_1 = tmp0_0*w2;
2215                            const double tmp8_1 = A_10*w1;
2216                            const double tmp9_1 = A_01*w2;
2217                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2218                            EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2219                            EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2220                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2221                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2222                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2223                            EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2224                            EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2225                            EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2226                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2227                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2228                            EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2229                            EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2230                            EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2231                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2232                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2233                        }
2234                        ///////////////
2235                        // process B //
2236                        ///////////////
2237                        if (!B.isEmpty()) {
2238                            add_EM_S=true;
2239                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2240                            const double tmp2_1 = B_p[0]*w8;
2241                            const double tmp0_1 = B_p[0]*w6;
2242                            const double tmp3_1 = B_p[1]*w9;
2243                            const double tmp1_1 = B_p[1]*w7;
2244                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2245                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2246                            EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
2247                            EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;
2248                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2249                            EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;
2250                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;
2251                            EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;
2252                            EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;
2253                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
2254                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;
2255                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2256                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;
2257                            EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;
2258                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;
2259                            EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
2260                        }
2261                        ///////////////
2262                        // process C //
2263                        ///////////////
2264                        if (!C.isEmpty()) {
2265                            add_EM_S=true;
2266                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2267                            const double tmp1_1 = C_p[1]*w7;
2268                            const double tmp0_1 = C_p[0]*w8;
2269                            const double tmp3_1 = C_p[0]*w6;
2270                            const double tmp2_1 = C_p[1]*w9;
2271                            EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2272                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2273                            EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
2274                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
2275                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2276                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;
2277                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;
2278                            EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;
2279                            EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;
2280                            EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
2281                            EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
2282                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2283                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;
2284                            EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;
2285                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2286                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;
2287                        }
2288                        ///////////////
2289                        // process D //
2290                        ///////////////
2291                        if (!D.isEmpty()) {
2292                            add_EM_S=true;
2293                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2294                            const double tmp0_1 = D_p[0]*w10;
2295                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
2296                            EM_S[INDEX2(1,0,4)]+=tmp0_1;
2297                            EM_S[INDEX2(2,0,4)]+=tmp0_1;
2298                            EM_S[INDEX2(3,0,4)]+=tmp0_1;
2299                            EM_S[INDEX2(0,1,4)]+=tmp0_1;
2300                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
2301                            EM_S[INDEX2(2,1,4)]+=tmp0_1;
2302                            EM_S[INDEX2(3,1,4)]+=tmp0_1;
2303                            EM_S[INDEX2(0,2,4)]+=tmp0_1;
2304                            EM_S[INDEX2(1,2,4)]+=tmp0_1;
2305                            EM_S[INDEX2(2,2,4)]+=tmp0_1;
2306                            EM_S[INDEX2(3,2,4)]+=tmp0_1;
2307                            EM_S[INDEX2(0,3,4)]+=tmp0_1;
2308                            EM_S[INDEX2(1,3,4)]+=tmp0_1;
2309                            EM_S[INDEX2(2,3,4)]+=tmp0_1;
2310                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
2311                        }
2312                        ///////////////
2313                        // process X //
2314                        ///////////////
2315                        if (!X.isEmpty()) {
2316                            add_EM_F=true;
2317                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2318                            const double tmp0_1 = X_p[0]*w11;
2319                            const double tmp2_1 = X_p[0]*w13;
2320                            const double tmp1_1 = X_p[1]*w12;
2321                            const double tmp3_1 = X_p[1]*w14;
2322                            EM_F[0]+=tmp0_1 + tmp1_1;
2323                            EM_F[1]+=tmp1_1 + tmp2_1;
2324                            EM_F[2]+=tmp0_1 + tmp3_1;
2325                            EM_F[3]+=tmp2_1 + tmp3_1;
2326                        }
2327                        ///////////////
2328                        // process Y //
2329                        ///////////////
2330                        if (!Y.isEmpty()) {
2331                            add_EM_F=true;
2332                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2333                            const double tmp0_1 = Y_p[0]*w15;
2334                            EM_F[0]+=tmp0_1;
2335                            EM_F[1]+=tmp0_1;
2336                            EM_F[2]+=tmp0_1;
2337                            EM_F[3]+=tmp0_1;
2338                        }
2339    
2340                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2341                        const index_t firstNode=m_N0*k1+k0;
2342                        IndexVector rowIndex;
2343                        rowIndex.push_back(m_dofMap[firstNode]);
2344                        rowIndex.push_back(m_dofMap[firstNode+1]);
2345                        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2346                        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2347                        if (add_EM_F) {
2348                            double *F_p=rhs.getSampleDataRW(0);
2349                            for (index_t i=0; i<rowIndex.size(); i++) {
2350                                if (rowIndex[i]<getNumDOF()) {
2351                                    F_p[rowIndex[i]]+=EM_F[i];
2352                                }
2353                            }
2354                        }
2355                        if (add_EM_S) {
2356                            addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2357                        }
2358    
2359                    } // end k0 loop
2360                } // end k1 loop
2361            } // end of colouring
2362      } // end of parallel region      } // end of parallel region
2363  }  }
2364    
# Line 2181  void Rectangle::assemblePDESingle(Paso_S Line 2366  void Rectangle::assemblePDESingle(Paso_S
2366  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2367          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2368          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2369          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  
2370  {  {
2371      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
2372      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 2194  void Rectangle::assemblePDESystem(Paso_S Line 2378  void Rectangle::assemblePDESystem(Paso_S
2378          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
2379      }      }
2380    
     /* GENERATOR SNIP_PDE_SYSTEM_PRE TOP */  
2381      const double w0 = -0.1555021169820365539*h1/h0;      const double w0 = -0.1555021169820365539*h1/h0;
2382      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
2383        const double w2 = -0.15550211698203655390;
2384        const double w3 = 0.041666666666666666667*h0/h1;
2385        const double w4 = 0.15550211698203655390;
2386        const double w5 = -0.041666666666666666667;
2387        const double w6 = -0.01116454968463011277*h1/h0;
2388        const double w7 = 0.011164549684630112770;
2389        const double w8 = -0.011164549684630112770;
2390        const double w9 = -0.041666666666666666667*h1/h0;
2391      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.041666666666666666667*h0/h1;
2392      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = 0.1555021169820365539*h1/h0;
2393      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = 0.1555021169820365539*h0/h1;
# Line 2207  void Rectangle::assemblePDESystem(Paso_S Line 2398  void Rectangle::assemblePDESystem(Paso_S
2398      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*h0/h1;
2399      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*h1/h0;
2400      const double w19 = 0.25000000000000000000;      const double w19 = 0.25000000000000000000;
     const double w2 = -0.15550211698203655390;  
2401      const double w20 = -0.25000000000000000000;      const double w20 = -0.25000000000000000000;
2402      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*h0/h1;
2403      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*h1/h0;
# Line 2218  void Rectangle::assemblePDESystem(Paso_S Line 2408  void Rectangle::assemblePDESystem(Paso_S
2408      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = -0.33333333333333333333*h0/h1;
2409      const double w28 = -0.032861463941450536761*h1;      const double w28 = -0.032861463941450536761*h1;
2410      const double w29 = -0.032861463941450536761*h0;      const double w29 = -0.032861463941450536761*h0;
     const double w3 = 0.041666666666666666667*h0/h1;  
2411      const double w30 = -0.12264065304058601714*h1;      const double w30 = -0.12264065304058601714*h1;
2412      const double w31 = -0.0023593469594139828636*h1;      const double w31 = -0.0023593469594139828636*h1;
2413      const double w32 = -0.008805202725216129906*h0;      const double w32 = -0.008805202725216129906*h0;
# Line 2229  void Rectangle::assemblePDESystem(Paso_S Line 2418  void Rectangle::assemblePDESystem(Paso_S
2418      const double w37 = 0.0023593469594139828636*h1;      const double w37 = 0.0023593469594139828636*h1;
2419      const double w38 = 0.12264065304058601714*h1;      const double w38 = 0.12264065304058601714*h1;
2420      const double w39 = 0.032861463941450536761*h0;      const double w39 = 0.032861463941450536761*h0;
     const double w4 = 0.15550211698203655390;  
2421      const double w40 = -0.12264065304058601714*h0;      const double w40 = -0.12264065304058601714*h0;
2422      const double w41 = -0.0023593469594139828636*h0;      const double w41 = -0.0023593469594139828636*h0;
2423      const double w42 = 0.0023593469594139828636*h0;      const double w42 = 0.0023593469594139828636*h0;
# Line 2240  void Rectangle::assemblePDESystem(Paso_S Line 2428  void Rectangle::assemblePDESystem(Paso_S
2428      const double w47 = 0.16666666666666666667*h1;      const double w47 = 0.16666666666666666667*h1;
2429      const double w48 = 0.083333333333333333333*h0;      const double w48 = 0.083333333333333333333*h0;
2430      const double w49 = -0.16666666666666666667*h0;      const double w49 = -0.16666666666666666667*h0;
     const double w5 = -0.041666666666666666667;  
2431      const double w50 = 0.16666666666666666667*h0;      const double w50 = 0.16666666666666666667*h0;
2432      const double w51 = -0.083333333333333333333*h1;      const double w51 = -0.083333333333333333333*h1;
2433      const double w52 = 0.025917019497006092316*h0*h1;      const double w52 = 0.025917019497006092316*h0*h1;
# Line 2251  void Rectangle::assemblePDESystem(Paso_S Line 2438  void Rectangle::assemblePDESystem(Paso_S
2438      const double w57 = 0.055555555555555555556*h0*h1;      const double w57 = 0.055555555555555555556*h0*h1;
2439      const double w58 = 0.027777777777777777778*h0*h1;      const double w58 = 0.027777777777777777778*h0*h1;
2440      const double w59 = 0.11111111111111111111*h0*h1;      const double w59 = 0.11111111111111111111*h0*h1;
     const double w6 = -0.01116454968463011277*h1/h0;  
2441      const double w60 = -0.19716878364870322056*h1;      const double w60 = -0.19716878364870322056*h1;
2442      const double w61 = -0.19716878364870322056*h0;      const double w61 = -0.19716878364870322056*h0;
2443      const double w62 = -0.052831216351296779436*h0;      const double w62 = -0.052831216351296779436*h0;
# Line 2262  void Rectangle::assemblePDESystem(Paso_S Line 2448  void Rectangle::assemblePDESystem(Paso_S
2448      const double w67 = 0.052831216351296779436*h0;      const double w67 = 0.052831216351296779436*h0;
2449      const double w68 = -0.5*h1;      const double w68 = -0.5*h1;
2450      const double w69 = -0.5*h0;      const double w69 = -0.5*h0;
     const double w7 = 0.011164549684630112770;  
2451      const double w70 = 0.5*h1;      const double w70 = 0.5*h1;
2452      const double w71 = 0.5*h0;      const double w71 = 0.5*h0;
2453      const double w72 = 0.1555021169820365539*h0*h1;      const double w72 = 0.1555021169820365539*h0*h1;
2454      const double w73 = 0.041666666666666666667*h0*h1;      const double w73 = 0.041666666666666666667*h0*h1;
2455      const double w74 = 0.01116454968463011277*h0*h1;      const double w74 = 0.01116454968463011277*h0*h1;
2456      const double w75 = 0.25*h0*h1;      const double w75 = 0.25*h0*h1;
     const double w8 = -0.011164549684630112770;  
     const double w9 = -0.041666666666666666667*h1/h0;  
     /* GENERATOR SNIP_PDE_SYSTEM_PRE BOTTOM */  
2457    
2458      rhs.requireWrite();      rhs.requireWrite();
2459  #pragma omp parallel  #pragma omp parallel
# Line 2285  void Rectangle::assemblePDESystem(Paso_S Line 2467  void Rectangle::assemblePDESystem(Paso_S
2467                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
2468                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
2469                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
                     /* GENERATOR SNIP_PDE_SYSTEM TOP */  
2470                      ///////////////                      ///////////////
2471                      // process A //                      // process A //
2472                      ///////////////                      ///////////////
# Line 2295  void Rectangle::assemblePDESystem(Paso_S Line 2476  void Rectangle::assemblePDESystem(Paso_S
2476                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
2477                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2478                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2479                                      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)];
2480                                      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)];
2481                                      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)];
2482                                      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)];
2483                                      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)];
2484                                      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)];
2485                                      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)];
2486                                      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)];
2487                                      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)];
2488                                      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)];
2489                                      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)];
2490                                      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)];
2491                                      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)];
2492                                      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)];
2493                                      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)];
2494                                      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)];
2495                                      const register double tmp4_0 = A_10_1 + A_10_2;                                      const double tmp0_0 = A_01_0 + A_01_3;
2496                                      const register double tmp12_0 = A_11_0 + A_11_2;                                      const double tmp1_0 = A_00_0 + A_00_1;
2497                                      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;
2498                                      const register double tmp10_0 = A_01_3 + A_10_3;                                      const double tmp3_0 = A_00_2 + A_00_3;
2499                                      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;
2500                                      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;
2501                                      const register double tmp13_0 = A_01_2 + A_10_1;                                      const double tmp6_0 = A_01_3 + A_10_0;
2502                                      const register double tmp3_0 = A_00_2 + A_00_3;                                      const double tmp7_0 = A_01_0 + A_10_3;
2503                                      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;
2504                                      const register double tmp18_0 = A_01_1 + A_10_1;                                      const double tmp9_0 = A_01_0 + A_10_0;
2505                                      const register double tmp1_0 = A_00_0 + A_00_1;                                      const double tmp10_0 = A_01_3 + A_10_3;
2506                                      const register double tmp15_0 = A_01_1 + A_10_2;                                      const double tmp11_0 = A_11_1 + A_11_3;
2507                                      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;
2508                                      const register double tmp16_0 = A_10_0 + A_10_3;                                      const double tmp13_0 = A_01_2 + A_10_1;
2509                                      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;
2510                                      const register double tmp17_0 = A_01_1 + A_01_2;                                      const double tmp15_0 = A_01_1 + A_10_2;
2511                                      const register double tmp9_0 = A_01_0 + A_10_0;                                      const double tmp16_0 = A_10_0 + A_10_3;
2512                                      const register double tmp7_0 = A_01_0 + A_10_3;                                      const double tmp17_0 = A_01_1 + A_01_2;
2513                                      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;
2514                                      const register double tmp19_0 = A_01_2 + A_10_2;                                      const double tmp19_0 = A_01_2 + A_10_2;
2515                                      const register double tmp14_1 = A_10_0*w8;                                      const double tmp0_1 = A_10_3*w8;
2516                                      const register double tmp23_1 = tmp3_0*w14;                                      const double tmp1_1 = tmp0_0*w1;
2517                                      const register double tmp35_1 = A_01_0*w8;                                      const double tmp2_1 = A_01_1*w4;
2518                                      const register double tmp54_1 = tmp13_0*w8;                                      const double tmp3_1 = tmp1_0*w0;
2519                                      const register double tmp20_1 = tmp9_0*w4;                                      const double tmp4_1 = A_01_2*w7;
2520                                      const register double tmp25_1 = tmp12_0*w12;                                      const double tmp5_1 = tmp2_0*w3;
2521                                      const register double tmp2_1 = A_01_1*w4;                                      const double tmp6_1 = tmp3_0*w6;
2522                                      const register double tmp44_1 = tmp7_0*w7;                                      const double tmp7_1 = A_10_0*w2;
2523                                      const register double tmp26_1 = tmp10_0*w4;                                      const double tmp8_1 = tmp4_0*w5;
2524                                      const register double tmp52_1 = tmp18_0*w8;                                      const double tmp9_1 = tmp2_0*w10;
2525                                      const register double tmp48_1 = A_10_1*w7;                                      const double tmp10_1 = tmp5_0*w9;
2526                                      const register double tmp46_1 = A_01_3*w8;                                      const double tmp11_1 = tmp6_0*w7;
2527                                      const register double tmp50_1 = A_01_0*w2;                                      const double tmp12_1 = tmp7_0*w4;
2528                                      const register double tmp8_1 = tmp4_0*w5;                                      const double tmp13_1 = tmp8_0*w1;
2529                                      const register double tmp56_1 = tmp19_0*w8;                                      const double tmp14_1 = A_10_0*w8;
2530                                      const register double tmp9_1 = tmp2_0*w10;                                      const double tmp15_1 = A_01_2*w4;
2531                                      const register double tmp19_1 = A_10_3*w2;                                      const double tmp16_1 = tmp3_0*w0;
2532                                      const register double tmp47_1 = A_10_2*w4;                                      const double tmp17_1 = A_01_1*w7;
2533                                      const register double tmp16_1 = tmp3_0*w0;                                      const double tmp18_1 = tmp1_0*w6;
2534                                      const register double tmp18_1 = tmp1_0*w6;                                      const double tmp19_1 = A_10_3*w2;
2535                                      const register double tmp31_1 = tmp11_0*w12;                                      const double tmp20_1 = tmp9_0*w4;
2536                                      const register double tmp55_1 = tmp15_0*w2;                                      const double tmp21_1 = tmp1_0*w11;
2537                                      const register double tmp39_1 = A_10_2*w7;                                      const double tmp22_1 = tmp10_0*w7;
2538                                      const register double tmp11_1 = tmp6_0*w7;                                      const double tmp23_1 = tmp3_0*w14;
2539                                      const register double tmp40_1 = tmp11_0*w17;                                      const double tmp24_1 = tmp11_0*w13;
2540                                      const register double tmp34_1 = tmp15_0*w8;                                      const double tmp25_1 = tmp12_0*w12;
2541                                      const register double tmp33_1 = tmp14_0*w5;                                      const double tmp26_1 = tmp10_0*w4;
2542                                      const register double tmp24_1 = tmp11_0*w13;                                      const double tmp27_1 = tmp3_0*w11;
2543                                      const register double tmp3_1 = tmp1_0*w0;                                      const double tmp28_1 = tmp9_0*w7;
2544                                      const register double tmp5_1 = tmp2_0*w3;                                      const double tmp29_1 = tmp1_0*w14;
2545                                      const register double tmp43_1 = tmp17_0*w5;                                      const double tmp30_1 = tmp12_0*w13;
2546                                      const register double tmp15_1 = A_01_2*w4;                                      const double tmp31_1 = tmp11_0*w12;
2547                                      const register double tmp53_1 = tmp19_0*w2;                                      const double tmp32_1 = tmp13_0*w2;
2548                                      const register double tmp27_1 = tmp3_0*w11;                                      const double tmp33_1 = tmp14_0*w5;
2549                                      const register double tmp32_1 = tmp13_0*w2;                                      const double tmp34_1 = tmp15_0*w8;
2550                                      const register double tmp10_1 = tmp5_0*w9;                                      const double tmp35_1 = A_01_0*w8;
2551                                      const register double tmp37_1 = A_10_1*w4;                                      const double tmp36_1 = tmp16_0*w1;
2552                                      const register double tmp38_1 = tmp5_0*w15;                                      const double tmp37_1 = A_10_1*w4;
2553                                      const register double tmp17_1 = A_01_1*w7;                                      const double tmp38_1 = tmp5_0*w15;
2554                                      const register double tmp12_1 = tmp7_0*w4;                                      const double tmp39_1 = A_10_2*w7;
2555                                      const register double tmp22_1 = tmp10_0*w7;                                      const double tmp40_1 = tmp11_0*w17;
2556                                      const register double tmp57_1 = tmp18_0*w2;                                      const double tmp41_1 = A_01_3*w2;
2557                                      const register double tmp28_1 = tmp9_0*w7;                                      const double tmp42_1 = tmp12_0*w16;
2558                                      const register double tmp29_1 = tmp1_0*w14;                                      const double tmp43_1 = tmp17_0*w5;
2559                                      const register double tmp51_1 = tmp11_0*w16;                                      const double tmp44_1 = tmp7_0*w7;
2560                                      const register double tmp42_1 = tmp12_0*w16;                                      const double tmp45_1 = tmp6_0*w4;
2561                                      const register double tmp49_1 = tmp12_0*w17;                                      const double tmp46_1 = A_01_3*w8;
2562                                      const register double tmp21_1 = tmp1_0*w11;                                      const double tmp47_1 = A_10_2*w4;
2563                                      const register double tmp1_1 = tmp0_0*w1;                                      const double tmp48_1 = A_10_1*w7;
2564                                      const register double tmp45_1 = tmp6_0*w4;                                      const double tmp49_1 = tmp12_0*w17;
2565                                      const register double tmp7_1 = A_10_0*w2;                                      const double tmp50_1 = A_01_0*w2;
2566                                      const register double tmp6_1 = tmp3_0*w6;                                      const double tmp51_1 = tmp11_0*w16;
2567                                      const register double tmp13_1 = tmp8_0*w1;                                      const double tmp52_1 = tmp18_0*w8;
2568                                      const register double tmp36_1 = tmp16_0*w1;                                      const double tmp53_1 = tmp19_0*w2;
2569                                      const register double tmp41_1 = A_01_3*w2;                                      const double tmp54_1 = tmp13_0*w8;
2570                                      const register double tmp30_1 = tmp12_0*w13;                                      const double tmp55_1 = tmp15_0*w2;
2571                                      const register double tmp4_1 = A_01_2*w7;                                      const double tmp56_1 = tmp19_0*w8;
2572                                      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;  
2573                                      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,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
2574                                      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;
2575                                        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;
2576                                      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,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
2577                                      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;
2578                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
2579                                      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,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
2580                                        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;
2581                                      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,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2582                                      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;  
2583                                      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,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
2584                                      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;
2585                                      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,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
2586                                      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;
2587                                        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;
2588                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2589                                  }                                  }
2590                              }                              }
2591                          } else { /* constant data */                          } else { // constant data
2592                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2593                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2594                                      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)];
2595                                      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)];
2596                                      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)];
2597                                      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)];
2598                                      const register double tmp0_0 = A_01 + A_10;                                      const double tmp0_0 = A_01 + A_10;
2599                                      const register double tmp0_1 = A_00*w18;                                      const double tmp0_1 = A_00*w18;
2600                                      const register double tmp10_1 = A_01*w20;                                      const double tmp1_1 = A_01*w19;
2601                                      const register double tmp12_1 = A_00*w26;                                      const double tmp2_1 = A_10*w20;
2602                                      const register double tmp4_1 = A_00*w22;                                      const double tmp3_1 = A_11*w21;
2603                                      const register double tmp8_1 = A_00*w24;                                      const double tmp4_1 = A_00*w22;
2604                                      const register double tmp13_1 = A_10*w19;                                      const double tmp5_1 = tmp0_0*w19;
2605                                      const register double tmp9_1 = tmp0_0*w20;                                      const double tmp6_1 = A_11*w23;
2606                                      const register double tmp3_1 = A_11*w21;                                      const double tmp7_1 = A_11*w25;
2607                                      const register double tmp11_1 = A_11*w27;                                      const double tmp8_1 = A_00*w24;
2608                                      const register double tmp1_1 = A_01*w19;                                      const double tmp9_1 = tmp0_0*w20;
2609                                      const register double tmp6_1 = A_11*w23;                                      const double tmp10_1 = A_01*w20;
2610                                      const register double tmp7_1 = A_11*w25;                                      const double tmp11_1 = A_11*w27;
2611                                      const register double tmp2_1 = A_10*w20;                                      const double tmp12_1 = A_00*w26;
2612                                      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;  
2613                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2614                                      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;
2615                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
2616                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2617                                      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;
2618                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2619                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
2620                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
2621                                      EM_S[INDEX4(k,m,0,2,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;
2622                                      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;  
2623                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2624                                      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;
2625                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2626                                      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;
2627                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2628                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2629                                  }                                  }
2630                              }                              }
2631                          }                          }
# Line 2458  void Rectangle::assemblePDESystem(Paso_S Line 2639  void Rectangle::assemblePDESystem(Paso_S
2639                          if (B.actsExpanded()) {                          if (B.actsExpanded()) {
2640                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2641                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2642                                      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)];
2643                                      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)];
2644                                      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)];
2645                                      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)];
2646                                      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)];
2647                                      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)];
2648                                      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)];
2649                                      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)];
2650                                      const register double tmp3_0 = B_0_0 + B_0_2;                                      const double tmp0_0 = B_1_0 + B_1_1;
2651                                      const register double tmp1_0 = B_1_2 + B_1_3;                                      const double tmp1_0 = B_1_2 + B_1_3;
2652                                      const register double tmp2_0 = B_0_1 + B_0_3;                                      const double tmp2_0 = B_0_1 + B_0_3;
2653                                      const register double tmp0_0 = B_1_0 + B_1_1;                                      const double tmp3_0 = B_0_0 + B_0_2;
2654                                      const register double tmp63_1 = B_1_1*w42;                                      const double tmp63_1 = B_1_1*w42;
2655                                      const register double tmp79_1 = B_1_1*w40;                                      const double tmp79_1 = B_1_1*w40;
2656                                      const register double tmp37_1 = tmp3_0*w35;                                      const double tmp37_1 = tmp3_0*w35;
2657                                      const register double tmp8_1 = tmp0_0*w32;                                      const double tmp8_1 = tmp0_0*w32;
2658                                      const register double tmp71_1 = B_0_1*w34;                                      const double tmp71_1 = B_0_1*w34;
2659                                      const register double tmp19_1 = B_0_3*w31;                                      const double tmp19_1 = B_0_3*w31;
2660                                      const register double tmp15_1 = B_0_3*w34;                                      const double tmp15_1 = B_0_3*w34;
2661                                      const register double tmp9_1 = tmp3_0*w34;                                      const double tmp9_1 = tmp3_0*w34;
2662                                      const register double tmp35_1 = B_1_0*w36;                                      const double tmp35_1 = B_1_0*w36;
2663                                      const register double tmp66_1 = B_0_3*w28;                                      const double tmp66_1 = B_0_3*w28;
2664                                      const register double tmp28_1 = B_1_0*w42;                                      const double tmp28_1 = B_1_0*w42;
2665                                      const register double tmp22_1 = B_1_0*w40;                                      const double tmp22_1 = B_1_0*w40;
2666                                      const register double tmp16_1 = B_1_2*w29;                                      const double tmp16_1 = B_1_2*w29;
2667                                      const register double tmp6_1 = tmp2_0*w35;                                      const double tmp6_1 = tmp2_0*w35;
2668                                      const register double tmp55_1 = B_1_3*w40;                                      const double tmp55_1 = B_1_3*w40;
2669                                      const register double tmp50_1 = B_1_3*w42;                                      const double tmp50_1 = B_1_3*w42;
2670                                      const register double tmp7_1 = tmp1_0*w29;                                      const double tmp7_1 = tmp1_0*w29;
2671                                      const register double tmp1_1 = tmp1_0*w32;                                      const double tmp1_1 = tmp1_0*w32;
2672                                      const register double tmp57_1 = B_0_3*w30;                                      const double tmp57_1 = B_0_3*w30;
2673                                      const register double tmp18_1 = B_1_1*w32;                                      const double tmp18_1 = B_1_1*w32;
2674                                      const register double tmp53_1 = B_1_0*w41;                                      const double tmp53_1 = B_1_0*w41;
2675                                      const register double tmp61_1 = B_1_3*w36;                                      const double tmp61_1 = B_1_3*w36;
2676                                      const register double tmp27_1 = B_0_3*w38;                                      const double tmp27_1 = B_0_3*w38;
2677                                      const register double tmp64_1 = B_0_2*w30;                                      const double tmp64_1 = B_0_2*w30;
2678                                      const register double tmp76_1 = B_0_1*w38;                                      const double tmp76_1 = B_0_1*w38;
2679                                      const register double tmp39_1 = tmp2_0*w34;                                      const double tmp39_1 = tmp2_0*w34;
2680                                      const register double tmp62_1 = B_0_1*w31;                                      const double tmp62_1 = B_0_1*w31;
2681                                      const register double tmp56_1 = B_0_0*w31;                                      const double tmp56_1 = B_0_0*w31;
2682                                      const register double tmp49_1 = B_1_1*w36;                                      const double tmp49_1 = B_1_1*w36;
2683                                      const register double tmp2_1 = B_0_2*w31;                                      const double tmp2_1 = B_0_2*w31;
2684                                      const register double tmp23_1 = B_0_2*w33;                                      const double tmp23_1 = B_0_2*w33;
2685                                      const register double tmp38_1 = B_1_1*w43;                                      const double tmp38_1 = B_1_1*w43;
2686                                      const register double tmp74_1 = B_1_2*w41;                                      const double tmp74_1 = B_1_2*w41;
2687                                      const register double tmp43_1 = B_1_1*w41;                                      const double tmp43_1 = B_1_1*w41;
2688                                      const register double tmp58_1 = B_0_2*w28;                                      const double tmp58_1 = B_0_2*w28;
2689                                      const register double tmp67_1 = B_0_0*w33;                                      const double tmp67_1 = B_0_0*w33;
2690                                      const register double tmp33_1 = tmp0_0*w39;                                      const double tmp33_1 = tmp0_0*w39;
2691                                      const register double tmp4_1 = B_0_0*w28;                                      const double tmp4_1 = B_0_0*w28;
2692                                      const register double tmp20_1 = B_0_0*w30;                                      const double tmp20_1 = B_0_0*w30;
2693                                      const register double tmp13_1 = B_0_2*w38;                                      const double tmp13_1 = B_0_2*w38;
2694                                      const register double tmp65_1 = B_1_2*w43;                                      const double tmp65_1 = B_1_2*w43;
2695                                      const register double tmp0_1 = tmp0_0*w29;                                      const double tmp0_1 = tmp0_0*w29;
2696                                      const register double tmp41_1 = tmp3_0*w33;                                      const double tmp41_1 = tmp3_0*w33;
2697                                      const register double tmp73_1 = B_0_2*w37;                                      const double tmp73_1 = B_0_2*w37;
2698                                      const register double tmp69_1 = B_0_0*w38;                                      const double tmp69_1 = B_0_0*w38;
2699                                      const register double tmp48_1 = B_1_2*w39;                                      const double tmp48_1 = B_1_2*w39;
2700                                      const register double tmp59_1 = B_0_1*w33;                                      const double tmp59_1 = B_0_1*w33;
2701                                      const register double tmp17_1 = B_1_3*w41;                                      const double tmp17_1 = B_1_3*w41;
2702                                      const register double tmp5_1 = B_0_3*w33;                                      const double tmp5_1 = B_0_3*w33;
2703                                      const register double tmp3_1 = B_0_1*w30;                                      const double tmp3_1 = B_0_1*w30;
2704                                      const register double tmp21_1 = B_0_1*w28;                                      const double tmp21_1 = B_0_1*w28;
2705                                      const register double tmp42_1 = B_1_0*w29;                                      const double tmp42_1 = B_1_0*w29;
2706                                      const register double tmp54_1 = B_1_2*w32;                                      const double tmp54_1 = B_1_2*w32;
2707                                      const register double tmp60_1 = B_1_0*w39;                                      const double tmp60_1 = B_1_0*w39;
2708                                      const register double tmp32_1 = tmp1_0*w36;                                      const double tmp32_1 = tmp1_0*w36;
2709                                      const register double tmp10_1 = B_0_1*w37;                                      const double tmp10_1 = B_0_1*w37;
2710                                      const register double tmp14_1 = B_0_0*w35;                                      const double tmp14_1 = B_0_0*w35;
2711                                      const register double tmp29_1 = B_0_1*w35;                                      const double tmp29_1 = B_0_1*w35;
2712                                      const register double tmp26_1 = B_1_2*w36;                                      const double tmp26_1 = B_1_2*w36;
2713                                      const register double tmp30_1 = B_1_3*w43;                                      const double tmp30_1 = B_1_3*w43;
2714                                      const register double tmp70_1 = B_0_2*w35;                                      const double tmp70_1 = B_0_2*w35;
2715                                      const register double tmp34_1 = B_1_3*w39;                                      const double tmp34_1 = B_1_3*w39;
2716                                      const register double tmp51_1 = B_1_0*w43;                                      const double tmp51_1 = B_1_0*w43;
2717                                      const register double tmp31_1 = B_0_2*w34;                                      const double tmp31_1 = B_0_2*w34;
2718                                      const register double tmp45_1 = tmp3_0*w28;                                      const double tmp45_1 = tmp3_0*w28;
2719                                      const register double tmp11_1 = tmp1_0*w39;                                      const double tmp11_1 = tmp1_0*w39;
2720                                      const register double tmp52_1 = B_1_1*w29;                                      const double tmp52_1 = B_1_1*w29;
2721                                      const register double tmp44_1 = B_1_3*w32;                                      const double tmp44_1 = B_1_3*w32;
2722                                      const register double tmp25_1 = B_1_1*w39;                                      const double tmp25_1 = B_1_1*w39;
2723                                      const register double tmp47_1 = tmp2_0*w33;                                      const double tmp47_1 = tmp2_0*w33;
2724                                      const register double tmp72_1 = B_1_3*w29;                                      const double tmp72_1 = B_1_3*w29;
2725                                      const register double tmp40_1 = tmp2_0*w28;                                      const double tmp40_1 = tmp2_0*w28;
2726                                      const register double tmp46_1 = B_1_2*w40;                                      const double tmp46_1 = B_1_2*w40;
2727                                      const register double tmp36_1 = B_1_2*w42;                                      const double tmp36_1 = B_1_2*w42;
2728                                      const register double tmp24_1 = B_0_0*w37;                                      const double tmp24_1 = B_0_0*w37;
2729                                      const register double tmp77_1 = B_0_3*w35;                                      const double tmp77_1 = B_0_3*w35;
2730                                      const register double tmp68_1 = B_0_3*w37;                                      const double tmp68_1 = B_0_3*w37;
2731                                      const register double tmp78_1 = B_0_0*w34;                                      const double tmp78_1 = B_0_0*w34;
2732                                      const register double tmp12_1 = tmp0_0*w36;                                      const double tmp12_1 = tmp0_0*w36;
2733                                      const register double tmp75_1 = B_1_0*w32;                                      const double tmp75_1 = B_1_0*w32;
2734                                      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,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2735                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2736                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_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 2749  void Rectangle::assemblePDESystem(Paso_S
2749                                      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;                                      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;
2750                                  }                                  }
2751                              }                              }
2752                          } else { /* constant data */                          } else { // constant data
2753                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2754                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2755                                      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)];
2756                                      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)];
2757                                      const register double tmp6_1 = B_1*w50;                                      const double tmp6_1 = B_1*w50;
2758                                      const register double tmp1_1 = B_1*w45;                                      const double tmp1_1 = B_1*w45;
2759                                      const register double tmp5_1 = B_1*w49;                                      const double tmp5_1 = B_1*w49;
2760                                      const register double tmp4_1 = B_1*w48;                                      const double tmp4_1 = B_1*w48;
2761                                      const register double tmp0_1 = B_0*w44;                                      const double tmp0_1 = B_0*w44;
2762                                      const register double tmp2_1 = B_0*w46;                                      const double tmp2_1 = B_0*w46;
2763                                      const register double tmp7_1 = B_0*w51;                                      const double tmp7_1 = B_0*w51;
2764                                      const register double tmp3_1 = B_0*w47;                                      const double tmp3_1 = B_0*w47;
2765                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
2766                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
2767                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
# Line 2610  void Rectangle::assemblePDESystem(Paso_S Line 2791  void Rectangle::assemblePDESystem(Paso_S
2791                          if (C.actsExpanded()) {                          if (C.actsExpanded()) {
2792                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2793                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2794                                      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)];
2795                                      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)];
2796                                      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)];
2797                                      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)];
2798                                      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)];
2799                                      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)];
2800                                      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)];
2801                                      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)];
2802                                      const register double tmp2_0 = C_0_1 + C_0_3;                                      const double tmp2_0 = C_0_1 + C_0_3;
2803                                      const register double tmp1_0 = C_1_2 + C_1_3;                                      const double tmp1_0 = C_1_2 + C_1_3;
2804                                      const register double tmp3_0 = C_0_0 + C_0_2;                                      const double tmp3_0 = C_0_0 + C_0_2;
2805                                      const register double tmp0_0 = C_1_0 + C_1_1;                                      const double tmp0_0 = C_1_0 + C_1_1;
2806                                      const register double tmp64_1 = C_0_2*w30;                                      const double tmp64_1 = C_0_2*w30;
2807                                      const register double tmp14_1 = C_0_2*w28;                                      const double tmp14_1 = C_0_2*w28;
2808                                      const register double tmp19_1 = C_0_3*w31;                                      const double tmp19_1 = C_0_3*w31;
2809                                      const register double tmp22_1 = C_1_0*w40;                                      const double tmp22_1 = C_1_0*w40;
2810                                      const register double tmp37_1 = tmp3_0*w35;                                      const double tmp37_1 = tmp3_0*w35;
2811                                      const register double tmp29_1 = C_0_1*w35;                                      const double tmp29_1 = C_0_1*w35;
2812                                      const register double tmp73_1 = C_0_2*w37;                                      const double tmp73_1 = C_0_2*w37;
2813                                      const register double tmp74_1 = C_1_2*w41;                                      const double tmp74_1 = C_1_2*w41;
2814                                      const register double tmp52_1 = C_1_3*w39;                                      const double tmp52_1 = C_1_3*w39;
2815                                      const register double tmp25_1 = C_1_1*w39;                                      const double tmp25_1 = C_1_1*w39;
2816                                      const register double tmp62_1 = C_0_1*w31;                                      const double tmp62_1 = C_0_1*w31;
2817                                      const register double tmp79_1 = C_1_1*w40;                                      const double tmp79_1 = C_1_1*w40;
2818                                      const register double tmp43_1 = C_1_1*w36;                                      const double tmp43_1 = C_1_1*w36;
2819                                      const register double tmp27_1 = C_0_3*w38;                                      const double tmp27_1 = C_0_3*w38;
2820                                      const register double tmp28_1 = C_1_0*w42;                                      const double tmp28_1 = C_1_0*w42;
2821                                      const register double tmp63_1 = C_1_1*w42;                                      const double tmp63_1 = C_1_1*w42;
2822                                      const register double tmp59_1 = C_0_3*w34;                                      const double tmp59_1 = C_0_3*w34;
2823                                      const register double tmp72_1 = C_1_3*w29;                                      const double tmp72_1 = C_1_3*w29;
2824                                      const register double tmp40_1 = tmp2_0*w35;                                      const double tmp40_1 = tmp2_0*w35;
2825                                      const register double tmp13_1 = C_0_3*w30;                                      const double tmp13_1 = C_0_3*w30;
2826                                      const register double tmp51_1 = C_1_2*w40;                                      const double tmp51_1 = C_1_2*w40;
2827                                      const register double tmp54_1 = C_1_2*w42;                                      const double tmp54_1 = C_1_2*w42;
2828                                      const register double tmp12_1 = C_0_0*w31;                                      const double tmp12_1 = C_0_0*w31;
2829                                      const register double tmp2_1 = tmp1_0*w32;                                      const double tmp2_1 = tmp1_0*w32;
2830                                      const register double tmp68_1 = C_0_2*w31;                                      const double tmp68_1 = C_0_2*w31;
2831                                      const register double tmp75_1 = C_1_0*w32;                                      const double tmp75_1 = C_1_0*w32;
2832                                      const register double tmp49_1 = C_1_1*w41;                                      const double tmp49_1 = C_1_1*w41;
2833                                      const register double tmp4_1 = C_0_2*w35;                                      const double tmp4_1 = C_0_2*w35;
2834                                      const register double tmp66_1 = C_0_3*w28;                                      const double tmp66_1 = C_0_3*w28;
2835                                      const register double tmp56_1 = C_0_1*w37;                                      const double tmp56_1 = C_0_1*w37;
2836                                      const register double tmp5_1 = C_0_1*w34;                                      const double tmp5_1 = C_0_1*w34;
2837                                      const register double tmp38_1 = tmp2_0*w34;                                      const double tmp38_1 = tmp2_0*w34;
2838                                      const register double tmp76_1 = C_0_1*w38;                                      const double tmp76_1 = C_0_1*w38;
2839                                      const register double tmp21_1 = C_0_1*w28;                                      const double tmp21_1 = C_0_1*w28;
2840                                      const register double tmp69_1 = C_0_1*w30;                                      const double tmp69_1 = C_0_1*w30;
2841                                      const register double tmp53_1 = C_1_0*w36;                                      const double tmp53_1 = C_1_0*w36;
2842                                      const register double tmp42_1 = C_1_2*w39;                                      const double tmp42_1 = C_1_2*w39;
2843                                      const register double tmp32_1 = tmp1_0*w29;                                      const double tmp32_1 = tmp1_0*w29;
2844                                      const register double tmp45_1 = C_1_0*w43;                                      const double tmp45_1 = C_1_0*w43;
2845                                      const register double tmp33_1 = tmp0_0*w32;                                      const double tmp33_1 = tmp0_0*w32;
2846                                      const register double tmp35_1 = C_1_0*w41;                                      const double tmp35_1 = C_1_0*w41;
2847                                      const register double tmp26_1 = C_1_2*w36;                                      const double tmp26_1 = C_1_2*w36;
2848                                      const register double tmp67_1 = C_0_0*w33;                                      const double tmp67_1 = C_0_0*w33;
2849                                      const register double tmp31_1 = C_0_2*w34;                                      const double tmp31_1 = C_0_2*w34;
2850                                      const register double tmp20_1 = C_0_0*w30;                                      const double tmp20_1 = C_0_0*w30;
2851                                      const register double tmp70_1 = C_0_0*w28;                                      const double tmp70_1 = C_0_0*w28;
2852                                      const register double tmp8_1 = tmp0_0*w39;                                      const double tmp8_1 = tmp0_0*w39;
2853                                      const register double tmp30_1 = C_1_3*w43;                                      const double tmp30_1 = C_1_3*w43;
2854                                      const register double tmp0_1 = tmp0_0*w29;                                      const double tmp0_1 = tmp0_0*w29;
2855                                      const register double tmp17_1 = C_1_3*w41;                                      const double tmp17_1 = C_1_3*w41;
2856                                      const register double tmp58_1 = C_0_0*w35;                                      const double tmp58_1 = C_0_0*w35;
2857                                      const register double tmp9_1 = tmp3_0*w33;                                      const double tmp9_1 = tmp3_0*w33;
2858                                      const register double tmp61_1 = C_1_3*w36;                                      const double tmp61_1 = C_1_3*w36;
2859                                      const register double tmp41_1 = tmp3_0*w34;                                      const double tmp41_1 = tmp3_0*w34;
2860                                      const register double tmp50_1 = C_1_3*w32;                                      const double tmp50_1 = C_1_3*w32;
2861                                      const register double tmp18_1 = C_1_1*w32;                                      const double tmp18_1 = C_1_1*w32;
2862                                      const register double tmp6_1 = tmp1_0*w36;                                      const double tmp6_1 = tmp1_0*w36;
2863                                      const register double tmp3_1 = C_0_0*w38;                                      const double tmp3_1 = C_0_0*w38;
2864                                      const register double tmp34_1 = C_1_1*w29;                                      const double tmp34_1 = C_1_1*w29;
2865                                      const register double tmp77_1 = C_0_3*w35;                                      const double tmp77_1 = C_0_3*w35;
2866                                      const register double tmp65_1 = C_1_2*w43;                                      const double tmp65_1 = C_1_2*w43;
2867                                      const register double tmp71_1 = C_0_3*w33;                                      const double tmp71_1 = C_0_3*w33;
2868                                      const register double tmp55_1 = C_1_1*w43;                                      const double tmp55_1 = C_1_1*w43;
2869                                      const register double tmp46_1 = tmp3_0*w28;                                      const double tmp46_1 = tmp3_0*w28;
2870                                      const register double tmp24_1 = C_0_0*w37;                                      const double tmp24_1 = C_0_0*w37;
2871                                      const register double tmp10_1 = tmp1_0*w39;                                      const double tmp10_1 = tmp1_0*w39;
2872                                      const register double tmp48_1 = C_1_0*w29;                                      const double tmp48_1 = C_1_0*w29;
2873                                      const register double tmp15_1 = C_0_1*w33;                                      const double tmp15_1 = C_0_1*w33;
2874                                      const register double tmp36_1 = C_1_2*w32;                                      const double tmp36_1 = C_1_2*w32;
2875                                      const register double tmp60_1 = C_1_0*w39;                                      const double tmp60_1 = C_1_0*w39;
2876                                      const register double tmp47_1 = tmp2_0*w33;                                      const double tmp47_1 = tmp2_0*w33;
2877                                      const register double tmp16_1 = C_1_2*w29;                                      const double tmp16_1 = C_1_2*w29;
2878                                      const register double tmp1_1 = C_0_3*w37;                                      const double tmp1_1 = C_0_3*w37;
2879                                      const register double tmp7_1 = tmp2_0*w28;                                      const double tmp7_1 = tmp2_0*w28;
2880                                      const register double tmp39_1 = C_1_3*w40;                                      const double tmp39_1 = C_1_3*w40;
2881                                      const register double tmp44_1 = C_1_3*w42;                                      const double tmp44_1 = C_1_3*w42;
2882                                      const register double tmp57_1 = C_0_2*w38;                                      const double tmp57_1 = C_0_2*w38;
2883                                      const register double tmp78_1 = C_0_0*w34;                                      const double tmp78_1 = C_0_0*w34;
2884                                      const register double tmp11_1 = tmp0_0*w36;                                      const double tmp11_1 = tmp0_0*w36;
2885                                      const register double tmp23_1 = C_0_2*w33;                                      const double tmp23_1 = C_0_2*w33;
2886                                      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,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2887                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2888                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_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 2901  void Rectangle::assemblePDESystem(Paso_S
2901                                      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;                                      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;
2902                                  }                                  }
2903                              }                              }
2904                          } else { /* constant data */                          } else { // constant data
2905                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2906                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2907                                      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)];
2908                                      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)];
2909                                      const register double tmp1_1 = C_1*w45;                                      const double tmp1_1 = C_1*w45;
2910                                      const register double tmp3_1 = C_0*w51;                                      const double tmp3_1 = C_0*w51;
2911                                      const register double tmp4_1 = C_0*w44;                                      const double tmp4_1 = C_0*w44;
2912                                      const register double tmp7_1 = C_0*w46;                                      const double tmp7_1 = C_0*w46;
2913                                      const register double tmp5_1 = C_1*w49;                                      const double tmp5_1 = C_1*w49;
2914                                      const register double tmp2_1 = C_1*w48;                                      const double tmp2_1 = C_1*w48;
2915                                      const register double tmp0_1 = C_0*w47;                                      const double tmp0_1 = C_0*w47;
2916                                      const register double tmp6_1 = C_1*w50;                                      const double tmp6_1 = C_1*w50;