/[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 3766 by caltinay, Thu Jan 12 08:17:49 2012 UTC revision 3776 by caltinay, Thu Jan 19 03:48:35 2012 UTC
# Line 76  Rectangle::Rectangle(int n0, int n1, dou Line 76  Rectangle::Rectangle(int n0, int n1, dou
76    
77      populateSampleIds();      populateSampleIds();
78      createPattern();      createPattern();
   
     if (m_gNE0>2*m_NX && m_gNE1>2*m_NY && m_gNE0%2==0 && m_gNE1%2==0) {  
         m_coarseMesh.reset(new Rectangle(n0/2, n1/2, l0, l1, d0, d1));  
     }  
79  }  }
80    
81  Rectangle::~Rectangle()  Rectangle::~Rectangle()
# Line 252  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
252              return &m_nodeId[0];              return &m_nodeId[0];
         case ReducedNodes:  
             if (m_coarseMesh)  
                 return m_coarseMesh->borrowSampleReferenceIDs(Nodes);  
             break;  
253          case DegreesOfFreedom:          case DegreesOfFreedom:
254            case ReducedDegreesOfFreedom: // FIXME: reduced
255              return &m_dofId[0];              return &m_dofId[0];
         case ReducedDegreesOfFreedom:  
             if (m_coarseMesh)  
                 return m_coarseMesh->borrowSampleReferenceIDs(DegreesOfFreedom);  
             break;  
256          case Elements:          case Elements:
257          case ReducedElements:          case ReducedElements:
258              return &m_elementId[0];              return &m_elementId[0];
# Line 286  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
280              return (m_dofMap[id] < getNumDOF());              return (m_dofMap[id] < getNumDOF());
         case ReducedNodes:  
             if (m_coarseMesh)  
                 return m_coarseMesh->ownSample(Nodes, id);  
             break;  
281          case DegreesOfFreedom:          case DegreesOfFreedom:
282          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
283              return true;              return true;
# Line 302  bool Rectangle::ownSample(int fsType, in Line 289  bool Rectangle::ownSample(int fsType, in
289          case ReducedFaceElements:          case ReducedFaceElements:
290              {              {
291                  // determine which face the sample belongs to before                  // determine which face the sample belongs to before
292                  // checking ownership of face element's first node                  // 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 310  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 501  void Rectangle::setToSize(escript::Data& Line 488  void Rectangle::setToSize(escript::Data&
488  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
489                                                  bool reducedColOrder) const                                                  bool reducedColOrder) const
490  {  {
491        /* FIXME: reduced
492      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
493          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
494        */
495      return m_pattern;      return m_pattern;
496  }  }
497    
# Line 643  void Rectangle::assembleGradient(escript Line 631  void Rectangle::assembleGradient(escript
631  #pragma omp parallel for  #pragma omp parallel for
632          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
633              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
634                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
635                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
636                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
637                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
638                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
639                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
640                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
# Line 665  void Rectangle::assembleGradient(escript Line 653  void Rectangle::assembleGradient(escript
653  #pragma omp parallel for  #pragma omp parallel for
654          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
655              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
656                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
657                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
658                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
659                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
660                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
661                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
662                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
# Line 683  void Rectangle::assembleGradient(escript Line 671  void Rectangle::assembleGradient(escript
671              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
672  #pragma omp for nowait  #pragma omp for nowait
673                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
674                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
675                      const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
676                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
677                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
678                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
679                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
680                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
# Line 699  void Rectangle::assembleGradient(escript Line 687  void Rectangle::assembleGradient(escript
687              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
688  #pragma omp for nowait  #pragma omp for nowait
689                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
690                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
691                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
692                      const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
693                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
694                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
695                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
696                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
# Line 715  void Rectangle::assembleGradient(escript Line 703  void Rectangle::assembleGradient(escript
703              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
704  #pragma omp for nowait  #pragma omp for nowait
705                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
706                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
707                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
708                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
709                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
710                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
711                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
712                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
# Line 731  void Rectangle::assembleGradient(escript Line 719  void Rectangle::assembleGradient(escript
719              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
720  #pragma omp for nowait  #pragma omp for nowait
721                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
722                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
723                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
724                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
725                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
726                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
727                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
728                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
# Line 752  void Rectangle::assembleGradient(escript Line 740  void Rectangle::assembleGradient(escript
740              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
741  #pragma omp for nowait  #pragma omp for nowait
742                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
743                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
744                      const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
745                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
746                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
747                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
748                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
749                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
# Line 766  void Rectangle::assembleGradient(escript Line 754  void Rectangle::assembleGradient(escript
754              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
755  #pragma omp for nowait  #pragma omp for nowait
756                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
757                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
758                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
759                      const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
760                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
761                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
762                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
763                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
# Line 780  void Rectangle::assembleGradient(escript Line 768  void Rectangle::assembleGradient(escript
768              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
769  #pragma omp for nowait  #pragma omp for nowait
770                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
771                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
772                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
773                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
774                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
775                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
776                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
777                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
# Line 794  void Rectangle::assembleGradient(escript Line 782  void Rectangle::assembleGradient(escript
782              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
783  #pragma omp for nowait  #pragma omp for nowait
784                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
785                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
786                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
787                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
788                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
789                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
790                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
791                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
# Line 827  void Rectangle::assembleIntegrate(vector Line 815  void Rectangle::assembleIntegrate(vector
815                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
816                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
817                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
818                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
819                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
820                          const register double f2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
821                          const register double f3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
822                          int_local[i]+=(f0+f1+f2+f3)*w;                          int_local[i]+=(f0+f1+f2+f3)*w;
823                      }  /* end of component loop i */                      }  /* end of component loop i */
824                  } /* end of k0 loop */                  } /* end of k0 loop */
# Line 870  void Rectangle::assembleIntegrate(vector Line 858  void Rectangle::assembleIntegrate(vector
858                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
859                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
860                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
861                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
862                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
863                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
864                      }  /* end of component loop i */                      }  /* end of component loop i */
865                  } /* end of k1 loop */                  } /* end of k1 loop */
# Line 882  void Rectangle::assembleIntegrate(vector Line 870  void Rectangle::assembleIntegrate(vector
870                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
871                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
872                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
873                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
874                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
875                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
876                      }  /* end of component loop i */                      }  /* end of component loop i */
877                  } /* end of k1 loop */                  } /* end of k1 loop */
# Line 894  void Rectangle::assembleIntegrate(vector Line 882  void Rectangle::assembleIntegrate(vector
882                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
883                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
884                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
885                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
886                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
887                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
888                      }  /* end of component loop i */                      }  /* end of component loop i */
889                  } /* end of k0 loop */                  } /* end of k0 loop */
# Line 906  void Rectangle::assembleIntegrate(vector Line 894  void Rectangle::assembleIntegrate(vector
894                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
895                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
896                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
897                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
898                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
899                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
900                      }  /* end of component loop i */                      }  /* end of component loop i */
901                  } /* end of k0 loop */                  } /* end of k0 loop */
# Line 1285  void Rectangle::createPattern() Line 1273  void Rectangle::createPattern()
1273      Paso_Pattern_free(rowPattern);      Paso_Pattern_free(rowPattern);
1274  }  }
1275    
1276    //private
1277    void Rectangle::addToMatrixAndRHS(Paso_SystemMatrix* S, escript::Data& F,
1278             const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
1279             bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
1280    {
1281        IndexVector rowIndex;
1282        rowIndex.push_back(m_dofMap[firstNode]);
1283        rowIndex.push_back(m_dofMap[firstNode+1]);
1284        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
1285        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
1286        if (addF) {
1287            double *F_p=F.getSampleDataRW(0);
1288            for (index_t i=0; i<rowIndex.size(); i++) {
1289                if (rowIndex[i]<getNumDOF()) {
1290                    for (index_t eq=0; eq<nEq; eq++) {
1291                        F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
1292                    }
1293                }
1294            }
1295        }
1296        if (addS) {
1297            addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
1298        }
1299    }
1300    
1301  //protected  //protected
1302  void Rectangle::interpolateNodesOnElements(escript::Data& out,  void Rectangle::interpolateNodesOnElements(escript::Data& out,
1303                                          escript::Data& in, bool reduced) const                                          escript::Data& in, bool reduced) const
# Line 1296  void Rectangle::interpolateNodesOnElemen Line 1309  void Rectangle::interpolateNodesOnElemen
1309  #pragma omp parallel for  #pragma omp parallel for
1310          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1311              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
1312                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1313                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1314                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1315                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1316                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1317                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1318                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
# Line 1314  void Rectangle::interpolateNodesOnElemen Line 1327  void Rectangle::interpolateNodesOnElemen
1327  #pragma omp parallel for  #pragma omp parallel for
1328          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1329              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
1330                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1331                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1332                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1333                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1334                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1335                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1336                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);
# Line 1343  void Rectangle::interpolateNodesOnFaces( Line 1356  void Rectangle::interpolateNodesOnFaces(
1356              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1357  #pragma omp for nowait  #pragma omp for nowait
1358                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1359                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1360                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1361                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1362                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1363                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
# Line 1354  void Rectangle::interpolateNodesOnFaces( Line 1367  void Rectangle::interpolateNodesOnFaces(
1367              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1368  #pragma omp for nowait  #pragma omp for nowait
1369                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1370                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1371                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1372                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1373                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1374                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
# Line 1365  void Rectangle::interpolateNodesOnFaces( Line 1378  void Rectangle::interpolateNodesOnFaces(
1378              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1379  #pragma omp for nowait  #pragma omp for nowait
1380                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1381                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1382                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1383                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1384                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1385                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
# Line 1376  void Rectangle::interpolateNodesOnFaces( Line 1389  void Rectangle::interpolateNodesOnFaces(
1389              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1390  #pragma omp for nowait  #pragma omp for nowait
1391                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1392                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1393                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1394                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1395                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1396                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
# Line 1394  void Rectangle::interpolateNodesOnFaces( Line 1407  void Rectangle::interpolateNodesOnFaces(
1407              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1408  #pragma omp for nowait  #pragma omp for nowait
1409                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1410                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1411                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1412                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1413                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1414                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;
# Line 1406  void Rectangle::interpolateNodesOnFaces( Line 1419  void Rectangle::interpolateNodesOnFaces(
1419              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1420  #pragma omp for nowait  #pragma omp for nowait
1421                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1422                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1423                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1424                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1425                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1426                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;
# Line 1418  void Rectangle::interpolateNodesOnFaces( Line 1431  void Rectangle::interpolateNodesOnFaces(
1431              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1432  #pragma omp for nowait  #pragma omp for nowait
1433                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1434                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1435                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1436                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1437                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1438                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;
# Line 1430  void Rectangle::interpolateNodesOnFaces( Line 1443  void Rectangle::interpolateNodesOnFaces(
1443              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1444  #pragma omp for nowait  #pragma omp for nowait
1445                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1446                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1447                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1448                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1449                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1450                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;
# Line 1447  void Rectangle::interpolateNodesOnFaces( Line 1460  void Rectangle::interpolateNodesOnFaces(
1460  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1461          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1462          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
1463          const escript::Data& X, const escript::Data& Y,          const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
1464  {  {
1465      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1466      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 1471  void Rectangle::assemblePDESingle(Paso_S Line 1483  void Rectangle::assemblePDESingle(Paso_S
1483      const double w16 = -0.01116454968463011277*h0/h1;      const double w16 = -0.01116454968463011277*h0/h1;
1484      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*h0/h1;
1485      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*h1/h0;
1486      const double w19 = 0.25000000000000000000;      const double w19 = 0.25;
1487      const double w20 = -0.25000000000000000000;      const double w20 = -0.25;
1488      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*h0/h1;
1489      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*h1/h0;
1490      const double w23 = -0.16666666666666666667*h0/h1;      const double w23 = -0.16666666666666666667*h0/h1;
# Line 1548  void Rectangle::assemblePDESingle(Paso_S Line 1560  void Rectangle::assemblePDESingle(Paso_S
1560                          add_EM_S=true;                          add_EM_S=true;
1561                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1562                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
1563                              const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];                              const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1564                              const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];                              const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1565                              const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];                              const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1566                              const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];                              const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1567                              const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];                              const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1568                              const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];                              const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1569                              const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];                              const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1570                              const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];                              const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1571                              const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];                              const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1572                              const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];                              const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1573                              const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];                              const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1574                              const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];                              const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1575                              const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];                              const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1576                              const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];                              const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1577                              const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];                              const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1578                              const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];                              const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1579                              const register double tmp0_0 = A_01_0 + A_01_3;                              const double tmp0_0 = A_01_0 + A_01_3;
1580                              const register double tmp1_0 = A_00_0 + A_00_1;                              const double tmp1_0 = A_00_0 + A_00_1;
1581                              const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                              const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1582                              const register double tmp3_0 = A_00_2 + A_00_3;                              const double tmp3_0 = A_00_2 + A_00_3;
1583                              const register double tmp4_0 = A_10_1 + A_10_2;                              const double tmp4_0 = A_10_1 + A_10_2;
1584                              const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                              const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1585                              const register double tmp6_0 = A_01_3 + A_10_0;                              const double tmp6_0 = A_01_3 + A_10_0;
1586                              const register double tmp7_0 = A_01_0 + A_10_3;                              const double tmp7_0 = A_01_0 + A_10_3;
1587                              const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                              const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1588                              const register double tmp9_0 = A_01_0 + A_10_0;                              const double tmp9_0 = A_01_0 + A_10_0;
1589                              const register double tmp12_0 = A_11_0 + A_11_2;                              const double tmp12_0 = A_11_0 + A_11_2;
1590                              const register double tmp10_0 = A_01_3 + A_10_3;                              const double tmp10_0 = A_01_3 + A_10_3;
1591                              const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                              const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
1592                              const register double tmp13_0 = A_01_2 + A_10_1;                              const double tmp13_0 = A_01_2 + A_10_1;
1593                              const register double tmp11_0 = A_11_1 + A_11_3;                              const double tmp11_0 = A_11_1 + A_11_3;
1594                              const register double tmp18_0 = A_01_1 + A_10_1;                              const double tmp18_0 = A_01_1 + A_10_1;
1595                              const register double tmp15_0 = A_01_1 + A_10_2;                              const double tmp15_0 = A_01_1 + A_10_2;
1596                              const register double tmp16_0 = A_10_0 + A_10_3;                              const double tmp16_0 = A_10_0 + A_10_3;
1597                              const register double tmp17_0 = A_01_1 + A_01_2;                              const double tmp17_0 = A_01_1 + A_01_2;
1598                              const register double tmp19_0 = A_01_2 + A_10_2;                              const double tmp19_0 = A_01_2 + A_10_2;
1599                              const register double tmp0_1 = A_10_3*w8;                              const double tmp0_1 = A_10_3*w8;
1600                              const register double tmp1_1 = tmp0_0*w1;                              const double tmp1_1 = tmp0_0*w1;
1601                              const register double tmp2_1 = A_01_1*w4;                              const double tmp2_1 = A_01_1*w4;
1602                              const register double tmp3_1 = tmp1_0*w0;                              const double tmp3_1 = tmp1_0*w0;
1603                              const register double tmp4_1 = A_01_2*w7;                              const double tmp4_1 = A_01_2*w7;
1604                              const register double tmp5_1 = tmp2_0*w3;                              const double tmp5_1 = tmp2_0*w3;
1605                              const register double tmp6_1 = tmp3_0*w6;                              const double tmp6_1 = tmp3_0*w6;
1606                              const register double tmp7_1 = A_10_0*w2;                              const double tmp7_1 = A_10_0*w2;
1607                              const register double tmp8_1 = tmp4_0*w5;                              const double tmp8_1 = tmp4_0*w5;
1608                              const register double tmp9_1 = tmp2_0*w10;                              const double tmp9_1 = tmp2_0*w10;
1609                              const register double tmp14_1 = A_10_0*w8;                              const double tmp14_1 = A_10_0*w8;
1610                              const register double tmp23_1 = tmp3_0*w14;                              const double tmp23_1 = tmp3_0*w14;
1611                              const register double tmp35_1 = A_01_0*w8;                              const double tmp35_1 = A_01_0*w8;
1612                              const register double tmp54_1 = tmp13_0*w8;                              const double tmp54_1 = tmp13_0*w8;
1613                              const register double tmp20_1 = tmp9_0*w4;                              const double tmp20_1 = tmp9_0*w4;
1614                              const register double tmp25_1 = tmp12_0*w12;                              const double tmp25_1 = tmp12_0*w12;
1615                              const register double tmp44_1 = tmp7_0*w7;                              const double tmp44_1 = tmp7_0*w7;
1616                              const register double tmp26_1 = tmp10_0*w4;                              const double tmp26_1 = tmp10_0*w4;
1617                              const register double tmp52_1 = tmp18_0*w8;                              const double tmp52_1 = tmp18_0*w8;
1618                              const register double tmp48_1 = A_10_1*w7;                              const double tmp48_1 = A_10_1*w7;
1619                              const register double tmp46_1 = A_01_3*w8;                              const double tmp46_1 = A_01_3*w8;
1620                              const register double tmp50_1 = A_01_0*w2;                              const double tmp50_1 = A_01_0*w2;
1621                              const register double tmp56_1 = tmp19_0*w8;                              const double tmp56_1 = tmp19_0*w8;
1622                              const register double tmp19_1 = A_10_3*w2;                              const double tmp19_1 = A_10_3*w2;
1623                              const register double tmp47_1 = A_10_2*w4;                              const double tmp47_1 = A_10_2*w4;
1624                              const register double tmp16_1 = tmp3_0*w0;                              const double tmp16_1 = tmp3_0*w0;
1625                              const register double tmp18_1 = tmp1_0*w6;                              const double tmp18_1 = tmp1_0*w6;
1626                              const register double tmp31_1 = tmp11_0*w12;                              const double tmp31_1 = tmp11_0*w12;
1627                              const register double tmp55_1 = tmp15_0*w2;                              const double tmp55_1 = tmp15_0*w2;
1628                              const register double tmp39_1 = A_10_2*w7;                              const double tmp39_1 = A_10_2*w7;
1629                              const register double tmp11_1 = tmp6_0*w7;                              const double tmp11_1 = tmp6_0*w7;
1630                              const register double tmp40_1 = tmp11_0*w17;                              const double tmp40_1 = tmp11_0*w17;
1631                              const register double tmp34_1 = tmp15_0*w8;                              const double tmp34_1 = tmp15_0*w8;
1632                              const register double tmp33_1 = tmp14_0*w5;                              const double tmp33_1 = tmp14_0*w5;
1633                              const register double tmp24_1 = tmp11_0*w13;                              const double tmp24_1 = tmp11_0*w13;
1634                              const register double tmp43_1 = tmp17_0*w5;                              const double tmp43_1 = tmp17_0*w5;
1635                              const register double tmp15_1 = A_01_2*w4;                              const double tmp15_1 = A_01_2*w4;
1636                              const register double tmp53_1 = tmp19_0*w2;                              const double tmp53_1 = tmp19_0*w2;
1637                              const register double tmp27_1 = tmp3_0*w11;                              const double tmp27_1 = tmp3_0*w11;
1638                              const register double tmp32_1 = tmp13_0*w2;                              const double tmp32_1 = tmp13_0*w2;
1639                              const register double tmp10_1 = tmp5_0*w9;                              const double tmp10_1 = tmp5_0*w9;
1640                              const register double tmp37_1 = A_10_1*w4;                              const double tmp37_1 = A_10_1*w4;
1641                              const register double tmp38_1 = tmp5_0*w15;                              const double tmp38_1 = tmp5_0*w15;
1642                              const register double tmp17_1 = A_01_1*w7;                              const double tmp17_1 = A_01_1*w7;
1643                              const register double tmp12_1 = tmp7_0*w4;                              const double tmp12_1 = tmp7_0*w4;
1644                              const register double tmp22_1 = tmp10_0*w7;                              const double tmp22_1 = tmp10_0*w7;
1645                              const register double tmp57_1 = tmp18_0*w2;                              const double tmp57_1 = tmp18_0*w2;
1646                              const register double tmp28_1 = tmp9_0*w7;                              const double tmp28_1 = tmp9_0*w7;
1647                              const register double tmp29_1 = tmp1_0*w14;                              const double tmp29_1 = tmp1_0*w14;
1648                              const register double tmp51_1 = tmp11_0*w16;                              const double tmp51_1 = tmp11_0*w16;
1649                              const register double tmp42_1 = tmp12_0*w16;                              const double tmp42_1 = tmp12_0*w16;
1650                              const register double tmp49_1 = tmp12_0*w17;                              const double tmp49_1 = tmp12_0*w17;
1651                              const register double tmp21_1 = tmp1_0*w11;                              const double tmp21_1 = tmp1_0*w11;
1652                              const register double tmp45_1 = tmp6_0*w4;                              const double tmp45_1 = tmp6_0*w4;
1653                              const register double tmp13_1 = tmp8_0*w1;                              const double tmp13_1 = tmp8_0*w1;
1654                              const register double tmp36_1 = tmp16_0*w1;                              const double tmp36_1 = tmp16_0*w1;
1655                              const register double tmp41_1 = A_01_3*w2;                              const double tmp41_1 = A_01_3*w2;
1656                              const register double tmp30_1 = tmp12_0*w13;                              const double tmp30_1 = tmp12_0*w13;
1657                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1658                              EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;                              EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
1659                              EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
# Line 1659  void Rectangle::assemblePDESingle(Paso_S Line 1671  void Rectangle::assemblePDESingle(Paso_S
1671                              EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;                              EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
1672                              EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1673                          } else { // constant data                          } else { // constant data
1674                              const register double A_00 = A_p[INDEX2(0,0,2)];                              const double A_00 = A_p[INDEX2(0,0,2)];
1675                              const register double A_10 = A_p[INDEX2(1,0,2)];                              const double A_10 = A_p[INDEX2(1,0,2)];
1676                              const register double A_01 = A_p[INDEX2(0,1,2)];                              const double A_01 = A_p[INDEX2(0,1,2)];
1677                              const register double A_11 = A_p[INDEX2(1,1,2)];                              const double A_11 = A_p[INDEX2(1,1,2)];
1678                              const register double tmp0_0 = A_01 + A_10;                              const double tmp0_0 = A_01 + A_10;
1679                              const register double tmp0_1 = A_00*w18;                              const double tmp0_1 = A_00*w18;
1680                              const register double tmp1_1 = A_01*w19;                              const double tmp1_1 = A_01*w19;
1681                              const register double tmp2_1 = A_10*w20;                              const double tmp2_1 = A_10*w20;
1682                              const register double tmp3_1 = A_11*w21;                              const double tmp3_1 = A_11*w21;
1683                              const register double tmp4_1 = A_00*w22;                              const double tmp4_1 = A_00*w22;
1684                              const register double tmp5_1 = tmp0_0*w19;                              const double tmp5_1 = tmp0_0*w19;
1685                              const register double tmp6_1 = A_11*w23;                              const double tmp6_1 = A_11*w23;
1686                              const register double tmp7_1 = A_11*w25;                              const double tmp7_1 = A_11*w25;
1687                              const register double tmp8_1 = A_00*w24;                              const double tmp8_1 = A_00*w24;
1688                              const register double tmp9_1 = tmp0_0*w20;                              const double tmp9_1 = tmp0_0*w20;
1689                              const register double tmp10_1 = A_01*w20;                              const double tmp10_1 = A_01*w20;
1690                              const register double tmp11_1 = A_11*w27;                              const double tmp11_1 = A_11*w27;
1691                              const register double tmp12_1 = A_00*w26;                              const double tmp12_1 = A_00*w26;
1692                              const register double tmp13_1 = A_10*w19;                              const double tmp13_1 = A_10*w19;
1693                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1694                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1695                              EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
# Line 1703  void Rectangle::assemblePDESingle(Paso_S Line 1715  void Rectangle::assemblePDESingle(Paso_S
1715                          add_EM_S=true;                          add_EM_S=true;
1716                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1717                          if (B.actsExpanded()) {                          if (B.actsExpanded()) {
1718                              const register double B_0_0 = B_p[INDEX2(0,0,2)];                              const double B_0_0 = B_p[INDEX2(0,0,2)];
1719                              const register double B_1_0 = B_p[INDEX2(1,0,2)];                              const double B_1_0 = B_p[INDEX2(1,0,2)];
1720                              const register double B_0_1 = B_p[INDEX2(0,1,2)];                              const double B_0_1 = B_p[INDEX2(0,1,2)];
1721                              const register double B_1_1 = B_p[INDEX2(1,1,2)];                              const double B_1_1 = B_p[INDEX2(1,1,2)];
1722                              const register double B_0_2 = B_p[INDEX2(0,2,2)];                              const double B_0_2 = B_p[INDEX2(0,2,2)];
1723                              const register double B_1_2 = B_p[INDEX2(1,2,2)];                              const double B_1_2 = B_p[INDEX2(1,2,2)];
1724                              const register double B_0_3 = B_p[INDEX2(0,3,2)];                              const double B_0_3 = B_p[INDEX2(0,3,2)];
1725                              const register double B_1_3 = B_p[INDEX2(1,3,2)];                              const double B_1_3 = B_p[INDEX2(1,3,2)];
1726                              const register double tmp0_0 = B_1_0 + B_1_1;                              const double tmp0_0 = B_1_0 + B_1_1;
1727                              const register double tmp1_0 = B_1_2 + B_1_3;                              const double tmp1_0 = B_1_2 + B_1_3;
1728                              const register double tmp2_0 = B_0_1 + B_0_3;                              const double tmp2_0 = B_0_1 + B_0_3;
1729                              const register double tmp3_0 = B_0_0 + B_0_2;                              const double tmp3_0 = B_0_0 + B_0_2;
1730                              const register double tmp63_1 = B_1_1*w42;                              const double tmp63_1 = B_1_1*w42;
1731                              const register double tmp79_1 = B_1_1*w40;                              const double tmp79_1 = B_1_1*w40;
1732                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
1733                              const register double tmp8_1 = tmp0_0*w32;                              const double tmp8_1 = tmp0_0*w32;
1734                              const register double tmp71_1 = B_0_1*w34;                              const double tmp71_1 = B_0_1*w34;
1735                              const register double tmp19_1 = B_0_3*w31;                              const double tmp19_1 = B_0_3*w31;
1736                              const register double tmp15_1 = B_0_3*w34;                              const double tmp15_1 = B_0_3*w34;
1737                              const register double tmp9_1 = tmp3_0*w34;                              const double tmp9_1 = tmp3_0*w34;
1738                              const register double tmp35_1 = B_1_0*w36;                              const double tmp35_1 = B_1_0*w36;
1739                              const register double tmp66_1 = B_0_3*w28;                              const double tmp66_1 = B_0_3*w28;
1740                              const register double tmp28_1 = B_1_0*w42;                              const double tmp28_1 = B_1_0*w42;
1741                              const register double tmp22_1 = B_1_0*w40;                              const double tmp22_1 = B_1_0*w40;
1742                              const register double tmp16_1 = B_1_2*w29;                              const double tmp16_1 = B_1_2*w29;
1743                              const register double tmp6_1 = tmp2_0*w35;                              const double tmp6_1 = tmp2_0*w35;
1744                              const register double tmp55_1 = B_1_3*w40;                              const double tmp55_1 = B_1_3*w40;
1745                              const register double tmp50_1 = B_1_3*w42;                              const double tmp50_1 = B_1_3*w42;
1746                              const register double tmp7_1 = tmp1_0*w29;                              const double tmp7_1 = tmp1_0*w29;
1747                              const register double tmp1_1 = tmp1_0*w32;                              const double tmp1_1 = tmp1_0*w32;
1748                              const register double tmp57_1 = B_0_3*w30;                              const double tmp57_1 = B_0_3*w30;
1749                              const register double tmp18_1 = B_1_1*w32;                              const double tmp18_1 = B_1_1*w32;
1750                              const register double tmp53_1 = B_1_0*w41;                              const double tmp53_1 = B_1_0*w41;
1751                              const register double tmp61_1 = B_1_3*w36;                              const double tmp61_1 = B_1_3*w36;
1752                              const register double tmp27_1 = B_0_3*w38;                              const double tmp27_1 = B_0_3*w38;
1753                              const register double tmp64_1 = B_0_2*w30;                              const double tmp64_1 = B_0_2*w30;
1754                              const register double tmp76_1 = B_0_1*w38;                              const double tmp76_1 = B_0_1*w38;
1755                              const register double tmp39_1 = tmp2_0*w34;                              const double tmp39_1 = tmp2_0*w34;
1756                              const register double tmp62_1 = B_0_1*w31;                              const double tmp62_1 = B_0_1*w31;
1757                              const register double tmp56_1 = B_0_0*w31;                              const double tmp56_1 = B_0_0*w31;
1758                              const register double tmp49_1 = B_1_1*w36;                              const double tmp49_1 = B_1_1*w36;
1759                              const register double tmp2_1 = B_0_2*w31;                              const double tmp2_1 = B_0_2*w31;
1760                              const register double tmp23_1 = B_0_2*w33;                              const double tmp23_1 = B_0_2*w33;
1761                              const register double tmp38_1 = B_1_1*w43;                              const double tmp38_1 = B_1_1*w43;
1762                              const register double tmp74_1 = B_1_2*w41;                              const double tmp74_1 = B_1_2*w41;
1763                              const register double tmp43_1 = B_1_1*w41;                              const double tmp43_1 = B_1_1*w41;
1764                              const register double tmp58_1 = B_0_2*w28;                              const double tmp58_1 = B_0_2*w28;
1765                              const register double tmp67_1 = B_0_0*w33;                              const double tmp67_1 = B_0_0*w33;
1766                              const register double tmp33_1 = tmp0_0*w39;                              const double tmp33_1 = tmp0_0*w39;
1767                              const register double tmp4_1 = B_0_0*w28;                              const double tmp4_1 = B_0_0*w28;
1768                              const register double tmp20_1 = B_0_0*w30;                              const double tmp20_1 = B_0_0*w30;
1769                              const register double tmp13_1 = B_0_2*w38;                              const double tmp13_1 = B_0_2*w38;
1770                              const register double tmp65_1 = B_1_2*w43;                              const double tmp65_1 = B_1_2*w43;
1771                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
1772                              const register double tmp41_1 = tmp3_0*w33;                              const double tmp41_1 = tmp3_0*w33;
1773                              const register double tmp73_1 = B_0_2*w37;                              const double tmp73_1 = B_0_2*w37;
1774                              const register double tmp69_1 = B_0_0*w38;                              const double tmp69_1 = B_0_0*w38;
1775                              const register double tmp48_1 = B_1_2*w39;                              const double tmp48_1 = B_1_2*w39;
1776                              const register double tmp59_1 = B_0_1*w33;                              const double tmp59_1 = B_0_1*w33;
1777                              const register double tmp17_1 = B_1_3*w41;                              const double tmp17_1 = B_1_3*w41;
1778                              const register double tmp5_1 = B_0_3*w33;                              const double tmp5_1 = B_0_3*w33;
1779                              const register double tmp3_1 = B_0_1*w30;                              const double tmp3_1 = B_0_1*w30;
1780                              const register double tmp21_1 = B_0_1*w28;                              const double tmp21_1 = B_0_1*w28;
1781                              const register double tmp42_1 = B_1_0*w29;                              const double tmp42_1 = B_1_0*w29;
1782                              const register double tmp54_1 = B_1_2*w32;                              const double tmp54_1 = B_1_2*w32;
1783                              const register double tmp60_1 = B_1_0*w39;                              const double tmp60_1 = B_1_0*w39;
1784                              const register double tmp32_1 = tmp1_0*w36;                              const double tmp32_1 = tmp1_0*w36;
1785                              const register double tmp10_1 = B_0_1*w37;                              const double tmp10_1 = B_0_1*w37;
1786                              const register double tmp14_1 = B_0_0*w35;                              const double tmp14_1 = B_0_0*w35;
1787                              const register double tmp29_1 = B_0_1*w35;                              const double tmp29_1 = B_0_1*w35;
1788                              const register double tmp26_1 = B_1_2*w36;                              const double tmp26_1 = B_1_2*w36;
1789                              const register double tmp30_1 = B_1_3*w43;                              const double tmp30_1 = B_1_3*w43;
1790                              const register double tmp70_1 = B_0_2*w35;                              const double tmp70_1 = B_0_2*w35;
1791                              const register double tmp34_1 = B_1_3*w39;                              const double tmp34_1 = B_1_3*w39;
1792                              const register double tmp51_1 = B_1_0*w43;                              const double tmp51_1 = B_1_0*w43;
1793                              const register double tmp31_1 = B_0_2*w34;                              const double tmp31_1 = B_0_2*w34;
1794                              const register double tmp45_1 = tmp3_0*w28;                              const double tmp45_1 = tmp3_0*w28;
1795                              const register double tmp11_1 = tmp1_0*w39;                              const double tmp11_1 = tmp1_0*w39;
1796                              const register double tmp52_1 = B_1_1*w29;                              const double tmp52_1 = B_1_1*w29;
1797                              const register double tmp44_1 = B_1_3*w32;                              const double tmp44_1 = B_1_3*w32;
1798                              const register double tmp25_1 = B_1_1*w39;                              const double tmp25_1 = B_1_1*w39;
1799                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
1800                              const register double tmp72_1 = B_1_3*w29;                              const double tmp72_1 = B_1_3*w29;
1801                              const register double tmp40_1 = tmp2_0*w28;                              const double tmp40_1 = tmp2_0*w28;
1802                              const register double tmp46_1 = B_1_2*w40;                              const double tmp46_1 = B_1_2*w40;
1803                              const register double tmp36_1 = B_1_2*w42;                              const double tmp36_1 = B_1_2*w42;
1804                              const register double tmp24_1 = B_0_0*w37;                              const double tmp24_1 = B_0_0*w37;
1805                              const register double tmp77_1 = B_0_3*w35;                              const double tmp77_1 = B_0_3*w35;
1806                              const register double tmp68_1 = B_0_3*w37;                              const double tmp68_1 = B_0_3*w37;
1807                              const register double tmp78_1 = B_0_0*w34;                              const double tmp78_1 = B_0_0*w34;
1808                              const register double tmp12_1 = tmp0_0*w36;                              const double tmp12_1 = tmp0_0*w36;
1809                              const register double tmp75_1 = B_1_0*w32;                              const double tmp75_1 = B_1_0*w32;
1810                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1811                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1812                              EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
# Line 1812  void Rectangle::assemblePDESingle(Paso_S Line 1824  void Rectangle::assemblePDESingle(Paso_S
1824                              EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;                              EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1825                              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(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1826                          } else { // constant data                          } else { // constant data
1827                              const register double B_0 = B_p[0];                              const double B_0 = B_p[0];
1828                              const register double B_1 = B_p[1];                              const double B_1 = B_p[1];
1829                              const register double tmp0_1 = B_0*w44;                              const double tmp0_1 = B_0*w44;
1830                              const register double tmp1_1 = B_1*w45;                              const double tmp1_1 = B_1*w45;
1831                              const register double tmp2_1 = B_0*w46;                              const double tmp2_1 = B_0*w46;
1832                              const register double tmp3_1 = B_0*w47;                              const double tmp3_1 = B_0*w47;
1833                              const register double tmp4_1 = B_1*w48;                              const double tmp4_1 = B_1*w48;
1834                              const register double tmp5_1 = B_1*w49;                              const double tmp5_1 = B_1*w49;
1835                              const register double tmp6_1 = B_1*w50;                              const double tmp6_1 = B_1*w50;
1836                              const register double tmp7_1 = B_0*w51;                              const double tmp7_1 = B_0*w51;
1837                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1838                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1839                              EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;                              EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
# Line 1847  void Rectangle::assemblePDESingle(Paso_S Line 1859  void Rectangle::assemblePDESingle(Paso_S
1859                          add_EM_S=true;                          add_EM_S=true;
1860                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
1861                          if (C.actsExpanded()) {                          if (C.actsExpanded()) {
1862                              const register double C_0_0 = C_p[INDEX2(0,0,2)];                              const double C_0_0 = C_p[INDEX2(0,0,2)];
1863                              const register double C_1_0 = C_p[INDEX2(1,0,2)];                              const double C_1_0 = C_p[INDEX2(1,0,2)];
1864                              const register double C_0_1 = C_p[INDEX2(0,1,2)];                              const double C_0_1 = C_p[INDEX2(0,1,2)];
1865                              const register double C_1_1 = C_p[INDEX2(1,1,2)];                              const double C_1_1 = C_p[INDEX2(1,1,2)];
1866                              const register double C_0_2 = C_p[INDEX2(0,2,2)];                              const double C_0_2 = C_p[INDEX2(0,2,2)];
1867                              const register double C_1_2 = C_p[INDEX2(1,2,2)];                              const double C_1_2 = C_p[INDEX2(1,2,2)];
1868                              const register double C_0_3 = C_p[INDEX2(0,3,2)];                              const double C_0_3 = C_p[INDEX2(0,3,2)];
1869                              const register double C_1_3 = C_p[INDEX2(1,3,2)];                              const double C_1_3 = C_p[INDEX2(1,3,2)];
1870                              const register double tmp0_0 = C_1_0 + C_1_1;                              const double tmp0_0 = C_1_0 + C_1_1;
1871                              const register double tmp1_0 = C_1_2 + C_1_3;                              const double tmp1_0 = C_1_2 + C_1_3;
1872                              const register double tmp2_0 = C_0_1 + C_0_3;                              const double tmp2_0 = C_0_1 + C_0_3;
1873                              const register double tmp3_0 = C_0_0 + C_0_2;                              const double tmp3_0 = C_0_0 + C_0_2;
1874                              const register double tmp64_1 = C_0_2*w30;                              const double tmp64_1 = C_0_2*w30;
1875                              const register double tmp14_1 = C_0_2*w28;                              const double tmp14_1 = C_0_2*w28;
1876                              const register double tmp19_1 = C_0_3*w31;                              const double tmp19_1 = C_0_3*w31;
1877                              const register double tmp22_1 = C_1_0*w40;                              const double tmp22_1 = C_1_0*w40;
1878                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
1879                              const register double tmp29_1 = C_0_1*w35;                              const double tmp29_1 = C_0_1*w35;
1880                              const register double tmp73_1 = C_0_2*w37;                              const double tmp73_1 = C_0_2*w37;
1881                              const register double tmp74_1 = C_1_2*w41;                              const double tmp74_1 = C_1_2*w41;
1882                              const register double tmp52_1 = C_1_3*w39;                              const double tmp52_1 = C_1_3*w39;
1883                              const register double tmp25_1 = C_1_1*w39;                              const double tmp25_1 = C_1_1*w39;
1884                              const register double tmp62_1 = C_0_1*w31;                              const double tmp62_1 = C_0_1*w31;
1885                              const register double tmp79_1 = C_1_1*w40;                              const double tmp79_1 = C_1_1*w40;
1886                              const register double tmp43_1 = C_1_1*w36;                              const double tmp43_1 = C_1_1*w36;
1887                              const register double tmp27_1 = C_0_3*w38;                              const double tmp27_1 = C_0_3*w38;
1888                              const register double tmp28_1 = C_1_0*w42;                              const double tmp28_1 = C_1_0*w42;
1889                              const register double tmp63_1 = C_1_1*w42;                              const double tmp63_1 = C_1_1*w42;
1890                              const register double tmp59_1 = C_0_3*w34;                              const double tmp59_1 = C_0_3*w34;
1891                              const register double tmp72_1 = C_1_3*w29;                              const double tmp72_1 = C_1_3*w29;
1892                              const register double tmp40_1 = tmp2_0*w35;                              const double tmp40_1 = tmp2_0*w35;
1893                              const register double tmp13_1 = C_0_3*w30;                              const double tmp13_1 = C_0_3*w30;
1894                              const register double tmp51_1 = C_1_2*w40;                              const double tmp51_1 = C_1_2*w40;
1895                              const register double tmp54_1 = C_1_2*w42;                              const double tmp54_1 = C_1_2*w42;
1896                              const register double tmp12_1 = C_0_0*w31;                              const double tmp12_1 = C_0_0*w31;
1897                              const register double tmp2_1 = tmp1_0*w32;                              const double tmp2_1 = tmp1_0*w32;
1898                              const register double tmp68_1 = C_0_2*w31;                              const double tmp68_1 = C_0_2*w31;
1899                              const register double tmp75_1 = C_1_0*w32;                              const double tmp75_1 = C_1_0*w32;
1900                              const register double tmp49_1 = C_1_1*w41;                              const double tmp49_1 = C_1_1*w41;
1901                              const register double tmp4_1 = C_0_2*w35;                              const double tmp4_1 = C_0_2*w35;
1902                              const register double tmp66_1 = C_0_3*w28;                              const double tmp66_1 = C_0_3*w28;
1903                              const register double tmp56_1 = C_0_1*w37;                              const double tmp56_1 = C_0_1*w37;
1904                              const register double tmp5_1 = C_0_1*w34;                              const double tmp5_1 = C_0_1*w34;
1905                              const register double tmp38_1 = tmp2_0*w34;                              const double tmp38_1 = tmp2_0*w34;
1906                              const register double tmp76_1 = C_0_1*w38;                              const double tmp76_1 = C_0_1*w38;
1907                              const register double tmp21_1 = C_0_1*w28;                              const double tmp21_1 = C_0_1*w28;
1908                              const register double tmp69_1 = C_0_1*w30;                              const double tmp69_1 = C_0_1*w30;
1909                              const register double tmp53_1 = C_1_0*w36;                              const double tmp53_1 = C_1_0*w36;
1910                              const register double tmp42_1 = C_1_2*w39;                              const double tmp42_1 = C_1_2*w39;
1911                              const register double tmp32_1 = tmp1_0*w29;                              const double tmp32_1 = tmp1_0*w29;
1912                              const register double tmp45_1 = C_1_0*w43;                              const double tmp45_1 = C_1_0*w43;
1913                              const register double tmp33_1 = tmp0_0*w32;                              const double tmp33_1 = tmp0_0*w32;
1914                              const register double tmp35_1 = C_1_0*w41;                              const double tmp35_1 = C_1_0*w41;
1915                              const register double tmp26_1 = C_1_2*w36;                              const double tmp26_1 = C_1_2*w36;
1916                              const register double tmp67_1 = C_0_0*w33;                              const double tmp67_1 = C_0_0*w33;
1917                              const register double tmp31_1 = C_0_2*w34;                              const double tmp31_1 = C_0_2*w34;
1918                              const register double tmp20_1 = C_0_0*w30;                              const double tmp20_1 = C_0_0*w30;
1919                              const register double tmp70_1 = C_0_0*w28;                              const double tmp70_1 = C_0_0*w28;
1920                              const register double tmp8_1 = tmp0_0*w39;                              const double tmp8_1 = tmp0_0*w39;
1921                              const register double tmp30_1 = C_1_3*w43;                              const double tmp30_1 = C_1_3*w43;
1922                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
1923                              const register double tmp17_1 = C_1_3*w41;                              const double tmp17_1 = C_1_3*w41;
1924                              const register double tmp58_1 = C_0_0*w35;                              const double tmp58_1 = C_0_0*w35;
1925                              const register double tmp9_1 = tmp3_0*w33;                              const double tmp9_1 = tmp3_0*w33;
1926                              const register double tmp61_1 = C_1_3*w36;                              const double tmp61_1 = C_1_3*w36;
1927                              const register double tmp41_1 = tmp3_0*w34;                              const double tmp41_1 = tmp3_0*w34;
1928                              const register double tmp50_1 = C_1_3*w32;                              const double tmp50_1 = C_1_3*w32;
1929                              const register double tmp18_1 = C_1_1*w32;                              const double tmp18_1 = C_1_1*w32;
1930                              const register double tmp6_1 = tmp1_0*w36;                              const double tmp6_1 = tmp1_0*w36;
1931                              const register double tmp3_1 = C_0_0*w38;                              const double tmp3_1 = C_0_0*w38;
1932                              const register double tmp34_1 = C_1_1*w29;                              const double tmp34_1 = C_1_1*w29;
1933                              const register double tmp77_1 = C_0_3*w35;                              const double tmp77_1 = C_0_3*w35;
1934                              const register double tmp65_1 = C_1_2*w43;                              const double tmp65_1 = C_1_2*w43;
1935                              const register double tmp71_1 = C_0_3*w33;                              const double tmp71_1 = C_0_3*w33;
1936                              const register double tmp55_1 = C_1_1*w43;                              const double tmp55_1 = C_1_1*w43;
1937                              const register double tmp46_1 = tmp3_0*w28;                              const double tmp46_1 = tmp3_0*w28;
1938                              const register double tmp24_1 = C_0_0*w37;                              const double tmp24_1 = C_0_0*w37;
1939                              const register double tmp10_1 = tmp1_0*w39;                              const double tmp10_1 = tmp1_0*w39;
1940                              const register double tmp48_1 = C_1_0*w29;                              const double tmp48_1 = C_1_0*w29;
1941                              const register double tmp15_1 = C_0_1*w33;                              const double tmp15_1 = C_0_1*w33;
1942                              const register double tmp36_1 = C_1_2*w32;                              const double tmp36_1 = C_1_2*w32;
1943                              const register double tmp60_1 = C_1_0*w39;                              const double tmp60_1 = C_1_0*w39;
1944                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
1945                              const register double tmp16_1 = C_1_2*w29;                              const double tmp16_1 = C_1_2*w29;
1946                              const register double tmp1_1 = C_0_3*w37;                              const double tmp1_1 = C_0_3*w37;
1947                              const register double tmp7_1 = tmp2_0*w28;                              const double tmp7_1 = tmp2_0*w28;
1948                              const register double tmp39_1 = C_1_3*w40;                              const double tmp39_1 = C_1_3*w40;
1949                              const register double tmp44_1 = C_1_3*w42;                              const double tmp44_1 = C_1_3*w42;
1950                              const register double tmp57_1 = C_0_2*w38;                              const double tmp57_1 = C_0_2*w38;
1951                              const register double tmp78_1 = C_0_0*w34;                              const double tmp78_1 = C_0_0*w34;
1952                              const register double tmp11_1 = tmp0_0*w36;                              const double tmp11_1 = tmp0_0*w36;
1953                              const register double tmp23_1 = C_0_2*w33;                              const double tmp23_1 = C_0_2*w33;
1954                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1955                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1956                              EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
# Line 1956  void Rectangle::assemblePDESingle(Paso_S Line 1968  void Rectangle::assemblePDESingle(Paso_S
1968                              EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;                              EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1969                              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(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1970                          } else { // constant data                          } else { // constant data
1971                              const register double C_0 = C_p[0];                              const double C_0 = C_p[0];
1972                              const register double C_1 = C_p[1];                              const double C_1 = C_p[1];
1973                              const register double tmp0_1 = C_0*w47;                              const double tmp0_1 = C_0*w47;
1974                              const register double tmp1_1 = C_1*w45;                              const double tmp1_1 = C_1*w45;
1975                              const register double tmp2_1 = C_1*w48;                              const double tmp2_1 = C_1*w48;
1976                              const register double tmp3_1 = C_0*w51;                              const double tmp3_1 = C_0*w51;
1977                              const register double tmp4_1 = C_0*w44;                              const double tmp4_1 = C_0*w44;
1978                              const register double tmp5_1 = C_1*w49;                              const double tmp5_1 = C_1*w49;
1979                              const register double tmp6_1 = C_1*w50;                              const double tmp6_1 = C_1*w50;
1980                              const register double tmp7_1 = C_0*w46;                              const double tmp7_1 = C_0*w46;
1981                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1982                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
1983                              EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
# Line 1991  void Rectangle::assemblePDESingle(Paso_S Line 2003  void Rectangle::assemblePDESingle(Paso_S
2003                          add_EM_S=true;                          add_EM_S=true;
2004                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2005                          if (D.actsExpanded()) {                          if (D.actsExpanded()) {
2006                              const register double D_0 = D_p[0];                              const double D_0 = D_p[0];
2007                              const register double D_1 = D_p[1];                              const double D_1 = D_p[1];
2008                              const register double D_2 = D_p[2];                              const double D_2 = D_p[2];
2009                              const register double D_3 = D_p[3];                              const double D_3 = D_p[3];
2010                              const register double tmp0_0 = D_0 + D_1;                              const double tmp0_0 = D_0 + D_1;
2011                              const register double tmp1_0 = D_2 + D_3;                              const double tmp1_0 = D_2 + D_3;
2012                              const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;                              const double tmp2_0 = D_0 + D_1 + D_2 + D_3;
2013                              const register double tmp3_0 = D_1 + D_2;                              const double tmp3_0 = D_1 + D_2;
2014                              const register double tmp4_0 = D_1 + D_3;                              const double tmp4_0 = D_1 + D_3;
2015                              const register double tmp5_0 = D_0 + D_2;                              const double tmp5_0 = D_0 + D_2;
2016                              const register double tmp6_0 = D_0 + D_3;                              const double tmp6_0 = D_0 + D_3;
2017                              const register double tmp0_1 = tmp0_0*w52;                              const double tmp0_1 = tmp0_0*w52;
2018                              const register double tmp1_1 = tmp1_0*w53;                              const double tmp1_1 = tmp1_0*w53;
2019                              const register double tmp2_1 = tmp2_0*w54;                              const double tmp2_1 = tmp2_0*w54;
2020                              const register double tmp3_1 = tmp1_0*w52;                              const double tmp3_1 = tmp1_0*w52;
2021                              const register double tmp4_1 = tmp0_0*w53;                              const double tmp4_1 = tmp0_0*w53;
2022                              const register double tmp5_1 = tmp3_0*w54;                              const double tmp5_1 = tmp3_0*w54;
2023                              const register double tmp6_1 = D_0*w55;                              const double tmp6_1 = D_0*w55;
2024                              const register double tmp7_1 = D_3*w56;                              const double tmp7_1 = D_3*w56;
2025                              const register double tmp8_1 = D_3*w55;                              const double tmp8_1 = D_3*w55;
2026                              const register double tmp9_1 = D_0*w56;                              const double tmp9_1 = D_0*w56;
2027                              const register double tmp10_1 = tmp4_0*w52;                              const double tmp10_1 = tmp4_0*w52;
2028                              const register double tmp11_1 = tmp5_0*w53;                              const double tmp11_1 = tmp5_0*w53;
2029                              const register double tmp12_1 = tmp5_0*w52;                              const double tmp12_1 = tmp5_0*w52;
2030                              const register double tmp13_1 = tmp4_0*w53;                              const double tmp13_1 = tmp4_0*w53;
2031                              const register double tmp14_1 = tmp6_0*w54;                              const double tmp14_1 = tmp6_0*w54;
2032                              const register double tmp15_1 = D_2*w55;                              const double tmp15_1 = D_2*w55;
2033                              const register double tmp16_1 = D_1*w56;                              const double tmp16_1 = D_1*w56;
2034                              const register double tmp17_1 = D_1*w55;                              const double tmp17_1 = D_1*w55;
2035                              const register double tmp18_1 = D_2*w56;                              const double tmp18_1 = D_2*w56;
2036                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2037                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2038                              EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;                              EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
# Line 2038  void Rectangle::assemblePDESingle(Paso_S Line 2050  void Rectangle::assemblePDESingle(Paso_S
2050                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2051                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2052                          } else { // constant data                          } else { // constant data
2053                              const register double tmp0_1 = D_p[0]*w57;                              const double tmp0_1 = D_p[0]*w57;
2054                              const register double tmp1_1 = D_p[0]*w58;                              const double tmp1_1 = D_p[0]*w58;
2055                              const register double tmp2_1 = D_p[0]*w59;                              const double tmp2_1 = D_p[0]*w59;
2056                              EM_S[INDEX2(0,0,4)]+=tmp2_1;                              EM_S[INDEX2(0,0,4)]+=tmp2_1;
2057                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
2058                              EM_S[INDEX2(2,0,4)]+=tmp0_1;                              EM_S[INDEX2(2,0,4)]+=tmp0_1;
# Line 2066  void Rectangle::assemblePDESingle(Paso_S Line 2078  void Rectangle::assemblePDESingle(Paso_S
2078                          add_EM_F=true;                          add_EM_F=true;
2079                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2080                          if (X.actsExpanded()) {                          if (X.actsExpanded()) {
2081                              const register double X_0_0 = X_p[INDEX2(0,0,2)];                              const double X_0_0 = X_p[INDEX2(0,0,2)];
2082                              const register double X_1_0 = X_p[INDEX2(1,0,2)];                              const double X_1_0 = X_p[INDEX2(1,0,2)];
2083                              const register double X_0_1 = X_p[INDEX2(0,1,2)];                              const double X_0_1 = X_p[INDEX2(0,1,2)];
2084                              const register double X_1_1 = X_p[INDEX2(1,1,2)];                              const double X_1_1 = X_p[INDEX2(1,1,2)];
2085                              const register double X_0_2 = X_p[INDEX2(0,2,2)];                              const double X_0_2 = X_p[INDEX2(0,2,2)];
2086                              const register double X_1_2 = X_p[INDEX2(1,2,2)];                              const double X_1_2 = X_p[INDEX2(1,2,2)];
2087                              const register double X_0_3 = X_p[INDEX2(0,3,2)];                              const double X_0_3 = X_p[INDEX2(0,3,2)];
2088                              const register double X_1_3 = X_p[INDEX2(1,3,2)];                              const double X_1_3 = X_p[INDEX2(1,3,2)];
2089                              const register double tmp0_0 = X_0_2 + X_0_3;                              const double tmp0_0 = X_0_2 + X_0_3;
2090                              const register double tmp1_0 = X_1_1 + X_1_3;                              const double tmp1_0 = X_1_1 + X_1_3;
2091                              const register double tmp2_0 = X_1_0 + X_1_2;                              const double tmp2_0 = X_1_0 + X_1_2;
2092                              const register double tmp3_0 = X_0_0 + X_0_1;                              const double tmp3_0 = X_0_0 + X_0_1;
2093                              const register double tmp0_1 = tmp0_0*w63;                              const double tmp0_1 = tmp0_0*w63;
2094                              const register double tmp1_1 = tmp1_0*w62;                              const double tmp1_1 = tmp1_0*w62;
2095                              const register double tmp2_1 = tmp2_0*w61;                              const double tmp2_1 = tmp2_0*w61;
2096                              const register double tmp3_1 = tmp3_0*w60;                              const double tmp3_1 = tmp3_0*w60;
2097                              const register double tmp4_1 = tmp0_0*w65;                              const double tmp4_1 = tmp0_0*w65;
2098                              const register double tmp5_1 = tmp3_0*w64;                              const double tmp5_1 = tmp3_0*w64;
2099                              const register double tmp6_1 = tmp2_0*w62;                              const double tmp6_1 = tmp2_0*w62;
2100                              const register double tmp7_1 = tmp1_0*w61;                              const double tmp7_1 = tmp1_0*w61;
2101                              const register double tmp8_1 = tmp2_0*w66;                              const double tmp8_1 = tmp2_0*w66;
2102                              const register double tmp9_1 = tmp3_0*w63;                              const double tmp9_1 = tmp3_0*w63;
2103                              const register double tmp10_1 = tmp0_0*w60;                              const double tmp10_1 = tmp0_0*w60;
2104                              const register double tmp11_1 = tmp1_0*w67;                              const double tmp11_1 = tmp1_0*w67;
2105                              const register double tmp12_1 = tmp1_0*w66;                              const double tmp12_1 = tmp1_0*w66;
2106                              const register double tmp13_1 = tmp3_0*w65;                              const double tmp13_1 = tmp3_0*w65;
2107                              const register double tmp14_1 = tmp0_0*w64;                              const double tmp14_1 = tmp0_0*w64;
2108                              const register double tmp15_1 = tmp2_0*w67;                              const double tmp15_1 = tmp2_0*w67;
2109                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2110                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2111                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2112                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2113                          } else { // constant data                          } else { // constant data
2114                              const register double tmp0_1 = X_p[1]*w69;                              const double tmp0_1 = X_p[1]*w69;
2115                              const register double tmp1_1 = X_p[0]*w68;                              const double tmp1_1 = X_p[0]*w68;
2116                              const register double tmp2_1 = X_p[0]*w70;                              const double tmp2_1 = X_p[0]*w70;
2117                              const register double tmp3_1 = X_p[1]*w71;                              const double tmp3_1 = X_p[1]*w71;
2118                              EM_F[0]+=tmp0_1 + tmp1_1;                              EM_F[0]+=tmp0_1 + tmp1_1;
2119                              EM_F[1]+=tmp0_1 + tmp2_1;                              EM_F[1]+=tmp0_1 + tmp2_1;
2120                              EM_F[2]+=tmp1_1 + tmp3_1;                              EM_F[2]+=tmp1_1 + tmp3_1;
# Line 2116  void Rectangle::assemblePDESingle(Paso_S Line 2128  void Rectangle::assemblePDESingle(Paso_S
2128                          add_EM_F=true;                          add_EM_F=true;
2129                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2130                          if (Y.actsExpanded()) {                          if (Y.actsExpanded()) {
2131                              const register double Y_0 = Y_p[0];                              const double Y_0 = Y_p[0];
2132                              const register double Y_1 = Y_p[1];                              const double Y_1 = Y_p[1];
2133                              const register double Y_2 = Y_p[2];                              const double Y_2 = Y_p[2];
2134                              const register double Y_3 = Y_p[3];                              const double Y_3 = Y_p[3];
2135                              const register double tmp0_0 = Y_1 + Y_2;                              const double tmp0_0 = Y_1 + Y_2;
2136                              const register double tmp1_0 = Y_0 + Y_3;                              const double tmp1_0 = Y_0 + Y_3;
2137                              const register double tmp0_1 = Y_0*w72;                              const double tmp0_1 = Y_0*w72;
2138                              const register double tmp1_1 = tmp0_0*w73;                              const double tmp1_1 = tmp0_0*w73;
2139                              const register double tmp2_1 = Y_3*w74;                              const double tmp2_1 = Y_3*w74;
2140                              const register double tmp3_1 = Y_1*w72;                              const double tmp3_1 = Y_1*w72;
2141                              const register double tmp4_1 = tmp1_0*w73;                              const double tmp4_1 = tmp1_0*w73;
2142                              const register double tmp5_1 = Y_2*w74;                              const double tmp5_1 = Y_2*w74;
2143                              const register double tmp6_1 = Y_2*w72;                              const double tmp6_1 = Y_2*w72;
2144                              const register double tmp7_1 = Y_1*w74;                              const double tmp7_1 = Y_1*w74;
2145                              const register double tmp8_1 = Y_3*w72;                              const double tmp8_1 = Y_3*w72;
2146                              const register double tmp9_1 = Y_0*w74;                              const double tmp9_1 = Y_0*w74;
2147                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2148                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2149                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2150                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2151                          } else { // constant data                          } else { // constant data
2152                              const register double tmp0_1 = Y_p[0]*w75;                              const double tmp0_1 = Y_p[0]*w75;
2153                              EM_F[0]+=tmp0_1;                              EM_F[0]+=tmp0_1;
2154                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
2155                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
# Line 2147  void Rectangle::assemblePDESingle(Paso_S Line 2159  void Rectangle::assemblePDESingle(Paso_S
2159    
2160                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2161                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
2162                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
                     rowIndex.push_back(m_dofMap[firstNode]);  
                     rowIndex.push_back(m_dofMap[firstNode+1]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                     if (add_EM_F) {  
                         double *F_p=rhs.getSampleDataRW(0);  
                         for (index_t i=0; i<rowIndex.size(); i++) {  
                             if (rowIndex[i]<getNumDOF()) {  
                                 F_p[rowIndex[i]]+=EM_F[i];  
                             }  
                         }  
                     }  
                     if (add_EM_S) {  
                         addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);  
                     }  
   
2163                  } // end k0 loop                  } // end k0 loop
2164              } // end k1 loop              } // end k1 loop
2165          } // end of colouring          } // end of colouring
# Line 2174  void Rectangle::assemblePDESingle(Paso_S Line 2170  void Rectangle::assemblePDESingle(Paso_S
2170  void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2171          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2172          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2173          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  
2174  {  {
2175      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
2176      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 2214  void Rectangle::assemblePDESingleReduced Line 2209  void Rectangle::assemblePDESingleReduced
2209                      if (!A.isEmpty()) {                      if (!A.isEmpty()) {
2210                          add_EM_S=true;                          add_EM_S=true;
2211                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2212                          const register double A_00 = A_p[INDEX2(0,0,2)];                          const double A_00 = A_p[INDEX2(0,0,2)];
2213                          const register double A_10 = A_p[INDEX2(1,0,2)];                          const double A_10 = A_p[INDEX2(1,0,2)];
2214                          const register double A_01 = A_p[INDEX2(0,1,2)];                          const double A_01 = A_p[INDEX2(0,1,2)];
2215                          const register double A_11 = A_p[INDEX2(1,1,2)];                          const double A_11 = A_p[INDEX2(1,1,2)];
2216                          const register double tmp0_0 = A_01 + A_10;                          const double tmp0_0 = A_01 + A_10;
2217                          const register double tmp0_1 = A_11*w3;                          const double tmp0_1 = A_11*w3;
2218                          const register double tmp1_1 = A_00*w0;                          const double tmp1_1 = A_00*w0;
2219                          const register double tmp2_1 = A_01*w1;                          const double tmp2_1 = A_01*w1;
2220                          const register double tmp3_1 = A_10*w2;                          const double tmp3_1 = A_10*w2;
2221                          const register double tmp4_1 = tmp0_0*w1;                          const double tmp4_1 = tmp0_0*w1;
2222                          const register double tmp5_1 = A_11*w4;                          const double tmp5_1 = A_11*w4;
2223                          const register double tmp6_1 = A_00*w5;                          const double tmp6_1 = A_00*w5;
2224                          const register double tmp7_1 = tmp0_0*w2;                          const double tmp7_1 = tmp0_0*w2;
2225                          const register double tmp8_1 = A_10*w1;                          const double tmp8_1 = A_10*w1;
2226                          const register double tmp9_1 = A_01*w2;                          const double tmp9_1 = A_01*w2;
2227                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2228                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2229                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
# Line 2252  void Rectangle::assemblePDESingleReduced Line 2247  void Rectangle::assemblePDESingleReduced
2247                      if (!B.isEmpty()) {                      if (!B.isEmpty()) {
2248                          add_EM_S=true;                          add_EM_S=true;
2249                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2250                          const register double tmp2_1 = B_p[0]*w8;                          const double tmp2_1 = B_p[0]*w8;
2251                          const register double tmp0_1 = B_p[0]*w6;                          const double tmp0_1 = B_p[0]*w6;
2252                          const register double tmp3_1 = B_p[1]*w9;                          const double tmp3_1 = B_p[1]*w9;
2253                          const register double tmp1_1 = B_p[1]*w7;                          const double tmp1_1 = B_p[1]*w7;
2254                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2255                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2256                          EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
# Line 2279  void Rectangle::assemblePDESingleReduced Line 2274  void Rectangle::assemblePDESingleReduced
2274                      if (!C.isEmpty()) {                      if (!C.isEmpty()) {
2275                          add_EM_S=true;                          add_EM_S=true;
2276                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2277                          const register double tmp1_1 = C_p[1]*w7;                          const double tmp1_1 = C_p[1]*w7;
2278                          const register double tmp0_1 = C_p[0]*w8;                          const double tmp0_1 = C_p[0]*w8;
2279                          const register double tmp3_1 = C_p[0]*w6;                          const double tmp3_1 = C_p[0]*w6;
2280                          const register double tmp2_1 = C_p[1]*w9;                          const double tmp2_1 = C_p[1]*w9;
2281                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2282                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2283                          EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
# Line 2306  void Rectangle::assemblePDESingleReduced Line 2301  void Rectangle::assemblePDESingleReduced
2301                      if (!D.isEmpty()) {                      if (!D.isEmpty()) {
2302                          add_EM_S=true;                          add_EM_S=true;
2303                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2304                          const register double tmp0_1 = D_p[0]*w10;                          const double tmp0_1 = D_p[0]*w10;
2305                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1;
2306                          EM_S[INDEX2(1,0,4)]+=tmp0_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1;
2307                          EM_S[INDEX2(2,0,4)]+=tmp0_1;                          EM_S[INDEX2(2,0,4)]+=tmp0_1;
# Line 2330  void Rectangle::assemblePDESingleReduced Line 2325  void Rectangle::assemblePDESingleReduced
2325                      if (!X.isEmpty()) {                      if (!X.isEmpty()) {
2326                          add_EM_F=true;                          add_EM_F=true;
2327                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2328                          const register double tmp0_1 = X_p[0]*w11;                          const double tmp0_1 = X_p[0]*w11;
2329                          const register double tmp2_1 = X_p[0]*w13;                          const double tmp2_1 = X_p[0]*w13;
2330                          const register double tmp1_1 = X_p[1]*w12;                          const double tmp1_1 = X_p[1]*w12;
2331                          const register double tmp3_1 = X_p[1]*w14;                          const double tmp3_1 = X_p[1]*w14;
2332                          EM_F[0]+=tmp0_1 + tmp1_1;                          EM_F[0]+=tmp0_1 + tmp1_1;
2333                          EM_F[1]+=tmp1_1 + tmp2_1;                          EM_F[1]+=tmp1_1 + tmp2_1;
2334                          EM_F[2]+=tmp0_1 + tmp3_1;                          EM_F[2]+=tmp0_1 + tmp3_1;
# Line 2345  void Rectangle::assemblePDESingleReduced Line 2340  void Rectangle::assemblePDESingleReduced
2340                      if (!Y.isEmpty()) {                      if (!Y.isEmpty()) {
2341                          add_EM_F=true;                          add_EM_F=true;
2342                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2343                          const register double tmp0_1 = Y_p[0]*w15;                          const double tmp0_1 = Y_p[0]*w15;
2344                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0_1;
2345                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0_1;
2346                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
# Line 2354  void Rectangle::assemblePDESingleReduced Line 2349  void Rectangle::assemblePDESingleReduced
2349    
2350                      // 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)
2351                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
2352                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
                     rowIndex.push_back(m_dofMap[firstNode]);  
                     rowIndex.push_back(m_dofMap[firstNode+1]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                     if (add_EM_F) {  
                         double *F_p=rhs.getSampleDataRW(0);  
                         for (index_t i=0; i<rowIndex.size(); i++) {  
                             if (rowIndex[i]<getNumDOF()) {  
                                 F_p[rowIndex[i]]+=EM_F[i];  
                             }  
                         }  
                     }  
                     if (add_EM_S) {  
                         addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);  
                     }  
   
2353                  } // end k0 loop                  } // end k0 loop
2354              } // end k1 loop              } // end k1 loop
2355          } // end of colouring          } // end of colouring
# Line 2381  void Rectangle::assemblePDESingleReduced Line 2360  void Rectangle::assemblePDESingleReduced
2360  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2361          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2362          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2363          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  
2364  {  {
2365      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
2366      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 2492  void Rectangle::assemblePDESystem(Paso_S Line 2470  void Rectangle::assemblePDESystem(Paso_S
2470                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
2471                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2472                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2473                                      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)];
2474                                      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)];
2475                                      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)];
2476                                      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)];
2477                                      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)];
2478                                      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)];
2479                                      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)];
2480                                      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)];
2481                                      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)];
2482                                      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)];
2483                                      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)];
2484                                      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)];
2485                                      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)];
2486                                      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)];
2487                                      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)];
2488                                      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)];
2489                                      const register double tmp0_0 = A_01_0 + A_01_3;                                      const double tmp0_0 = A_01_0 + A_01_3;
2490                                      const register double tmp1_0 = A_00_0 + A_00_1;                                      const double tmp1_0 = A_00_0 + A_00_1;
2491                                      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;
2492                                      const register double tmp3_0 = A_00_2 + A_00_3;                                      const double tmp3_0 = A_00_2 + A_00_3;
2493                                      const register double tmp4_0 = A_10_1 + A_10_2;                                      const double tmp4_0 = A_10_1 + A_10_2;
2494                                      const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                                      const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
2495                                      const register double tmp6_0 = A_01_3 + A_10_0;                                      const double tmp6_0 = A_01_3 + A_10_0;
2496                                      const register double tmp7_0 = A_01_0 + A_10_3;                                      const double tmp7_0 = A_01_0 + A_10_3;
2497                                      const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                                      const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
2498                                      const register double tmp9_0 = A_01_0 + A_10_0;                                      const double tmp9_0 = A_01_0 + A_10_0;
2499                                      const register double tmp10_0 = A_01_3 + A_10_3;                                      const double tmp10_0 = A_01_3 + A_10_3;
2500                                      const register double tmp11_0 = A_11_1 + A_11_3;                                      const double tmp11_0 = A_11_1 + A_11_3;
2501                                      const register double tmp12_0 = A_11_0 + A_11_2;                                      const double tmp12_0 = A_11_0 + A_11_2;
2502                                      const register double tmp13_0 = A_01_2 + A_10_1;                                      const double tmp13_0 = A_01_2 + A_10_1;
2503                                      const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                                      const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
2504                                      const register double tmp15_0 = A_01_1 + A_10_2;                                      const double tmp15_0 = A_01_1 + A_10_2;
2505                                      const register double tmp16_0 = A_10_0 + A_10_3;                                      const double tmp16_0 = A_10_0 + A_10_3;
2506                                      const register double tmp17_0 = A_01_1 + A_01_2;                                      const double tmp17_0 = A_01_1 + A_01_2;
2507                                      const register double tmp18_0 = A_01_1 + A_10_1;                                      const double tmp18_0 = A_01_1 + A_10_1;
2508                                      const register double tmp19_0 = A_01_2 + A_10_2;                                      const double tmp19_0 = A_01_2 + A_10_2;
2509                                      const register double tmp0_1 = A_10_3*w8;                                      const double tmp0_1 = A_10_3*w8;
2510                                      const register double tmp1_1 = tmp0_0*w1;                                      const double tmp1_1 = tmp0_0*w1;
2511                                      const register double tmp2_1 = A_01_1*w4;                                      const double tmp2_1 = A_01_1*w4;
2512                                      const register double tmp3_1 = tmp1_0*w0;                                      const double tmp3_1 = tmp1_0*w0;
2513                                      const register double tmp4_1 = A_01_2*w7;                                      const double tmp4_1 = A_01_2*w7;
2514                                      const register double tmp5_1 = tmp2_0*w3;                                      const double tmp5_1 = tmp2_0*w3;
2515                                      const register double tmp6_1 = tmp3_0*w6;                                      const double tmp6_1 = tmp3_0*w6;
2516                                      const register double tmp7_1 = A_10_0*w2;                                      const double tmp7_1 = A_10_0*w2;
2517                                      const register double tmp8_1 = tmp4_0*w5;                                      const double tmp8_1 = tmp4_0*w5;
2518                                      const register double tmp9_1 = tmp2_0*w10;                                      const double tmp9_1 = tmp2_0*w10;
2519                                      const register double tmp10_1 = tmp5_0*w9;                                      const double tmp10_1 = tmp5_0*w9;
2520                                      const register double tmp11_1 = tmp6_0*w7;                                      const double tmp11_1 = tmp6_0*w7;
2521                                      const register double tmp12_1 = tmp7_0*w4;                                      const double tmp12_1 = tmp7_0*w4;
2522                                      const register double tmp13_1 = tmp8_0*w1;                                      const double tmp13_1 = tmp8_0*w1;
2523                                      const register double tmp14_1 = A_10_0*w8;                                      const double tmp14_1 = A_10_0*w8;
2524                                      const register double tmp15_1 = A_01_2*w4;                                      const double tmp15_1 = A_01_2*w4;
2525                                      const register double tmp16_1 = tmp3_0*w0;                                      const double tmp16_1 = tmp3_0*w0;
2526                                      const register double tmp17_1 = A_01_1*w7;                                      const double tmp17_1 = A_01_1*w7;
2527                                      const register double tmp18_1 = tmp1_0*w6;                                      const double tmp18_1 = tmp1_0*w6;
2528                                      const register double tmp19_1 = A_10_3*w2;                                      const double tmp19_1 = A_10_3*w2;
2529                                      const register double tmp20_1 = tmp9_0*w4;                                      const double tmp20_1 = tmp9_0*w4;
2530                                      const register double tmp21_1 = tmp1_0*w11;                                      const double tmp21_1 = tmp1_0*w11;
2531                                      const register double tmp22_1 = tmp10_0*w7;                                      const double tmp22_1 = tmp10_0*w7;
2532                                      const register double tmp23_1 = tmp3_0*w14;                                      const double tmp23_1 = tmp3_0*w14;
2533                                      const register double tmp24_1 = tmp11_0*w13;                                      const double tmp24_1 = tmp11_0*w13;
2534                                      const register double tmp25_1 = tmp12_0*w12;                                      const double tmp25_1 = tmp12_0*w12;
2535                                      const register double tmp26_1 = tmp10_0*w4;                                      const double tmp26_1 = tmp10_0*w4;
2536                                      const register double tmp27_1 = tmp3_0*w11;                                      const double tmp27_1 = tmp3_0*w11;
2537                                      const register double tmp28_1 = tmp9_0*w7;                                      const double tmp28_1 = tmp9_0*w7;
2538                                      const register double tmp29_1 = tmp1_0*w14;                                      const double tmp29_1 = tmp1_0*w14;
2539                                      const register double tmp30_1 = tmp12_0*w13;                                      const double tmp30_1 = tmp12_0*w13;
2540                                      const register double tmp31_1 = tmp11_0*w12;                                      const double tmp31_1 = tmp11_0*w12;
2541                                      const register double tmp32_1 = tmp13_0*w2;                                      const double tmp32_1 = tmp13_0*w2;
2542                                      const register double tmp33_1 = tmp14_0*w5;                                      const double tmp33_1 = tmp14_0*w5;
2543                                      const register double tmp34_1 = tmp15_0*w8;                                      const double tmp34_1 = tmp15_0*w8;
2544                                      const register double tmp35_1 = A_01_0*w8;                                      const double tmp35_1 = A_01_0*w8;
2545                                      const register double tmp36_1 = tmp16_0*w1;                                      const double tmp36_1 = tmp16_0*w1;
2546                                      const register double tmp37_1 = A_10_1*w4;                                      const double tmp37_1 = A_10_1*w4;
2547                                      const register double tmp38_1 = tmp5_0*w15;                                      const double tmp38_1 = tmp5_0*w15;
2548                                      const register double tmp39_1 = A_10_2*w7;                                      const double tmp39_1 = A_10_2*w7;
2549                                      const register double tmp40_1 = tmp11_0*w17;                                      const double tmp40_1 = tmp11_0*w17;
2550                                      const register double tmp41_1 = A_01_3*w2;                                      const double tmp41_1 = A_01_3*w2;
2551                                      const register double tmp42_1 = tmp12_0*w16;                                      const double tmp42_1 = tmp12_0*w16;
2552                                      const register double tmp43_1 = tmp17_0*w5;                                      const double tmp43_1 = tmp17_0*w5;
2553                                      const register double tmp44_1 = tmp7_0*w7;                                      const double tmp44_1 = tmp7_0*w7;
2554                                      const register double tmp45_1 = tmp6_0*w4;                                      const double tmp45_1 = tmp6_0*w4;
2555                                      const register double tmp46_1 = A_01_3*w8;                                      const double tmp46_1 = A_01_3*w8;
2556                                      const register double tmp47_1 = A_10_2*w4;                                      const double tmp47_1 = A_10_2*w4;
2557                                      const register double tmp48_1 = A_10_1*w7;                                      const double tmp48_1 = A_10_1*w7;
2558                                      const register double tmp49_1 = tmp12_0*w17;                                      const double tmp49_1 = tmp12_0*w17;
2559                                      const register double tmp50_1 = A_01_0*w2;                                      const double tmp50_1 = A_01_0*w2;
2560                                      const register double tmp51_1 = tmp11_0*w16;                                      const double tmp51_1 = tmp11_0*w16;
2561                                      const register double tmp52_1 = tmp18_0*w8;                                      const double tmp52_1 = tmp18_0*w8;
2562                                      const register double tmp53_1 = tmp19_0*w2;                                      const double tmp53_1 = tmp19_0*w2;
2563                                      const register double tmp54_1 = tmp13_0*w8;                                      const double tmp54_1 = tmp13_0*w8;
2564                                      const register double tmp55_1 = tmp15_0*w2;                                      const double tmp55_1 = tmp15_0*w2;
2565                                      const register double tmp56_1 = tmp19_0*w8;                                      const double tmp56_1 = tmp19_0*w8;
2566                                      const register double tmp57_1 = tmp18_0*w2;                                      const double tmp57_1 = tmp18_0*w2;
2567                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
2568                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
2569                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;                                      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;
# Line 2607  void Rectangle::assemblePDESystem(Paso_S Line 2585  void Rectangle::assemblePDESystem(Paso_S
2585                          } else { // constant data                          } else { // constant data
2586                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2587                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2588                                      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)];
2589                                      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)];
2590                                      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)];
2591                                      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)];
2592                                      const register double tmp0_0 = A_01 + A_10;                                      const double tmp0_0 = A_01 + A_10;
2593                                      const register double tmp0_1 = A_00*w18;                                      const double tmp0_1 = A_00*w18;
2594                                      const register double tmp1_1 = A_01*w19;                                      const double tmp1_1 = A_01*w19;
2595                                      const register double tmp2_1 = A_10*w20;                                      const double tmp2_1 = A_10*w20;
2596                                      const register double tmp3_1 = A_11*w21;                                      const double tmp3_1 = A_11*w21;
2597                                      const register double tmp4_1 = A_00*w22;                                      const double tmp4_1 = A_00*w22;
2598                                      const register double tmp5_1 = tmp0_0*w19;                                      const double tmp5_1 = tmp0_0*w19;
2599                                      const register double tmp6_1 = A_11*w23;                                      const double tmp6_1 = A_11*w23;
2600                                      const register double tmp7_1 = A_11*w25;                                      const double tmp7_1 = A_11*w25;
2601                                      const register double tmp8_1 = A_00*w24;                                      const double tmp8_1 = A_00*w24;
2602                                      const register double tmp9_1 = tmp0_0*w20;                                      const double tmp9_1 = tmp0_0*w20;
2603                                      const register double tmp10_1 = A_01*w20;                                      const double tmp10_1 = A_01*w20;
2604                                      const register double tmp11_1 = A_11*w27;                                      const double tmp11_1 = A_11*w27;
2605                                      const register double tmp12_1 = A_00*w26;                                      const double tmp12_1 = A_00*w26;
2606                                      const register double tmp13_1 = A_10*w19;                                      const double tmp13_1 = A_10*w19;
2607                                      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;
2608                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2609                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
# Line 2655  void Rectangle::assemblePDESystem(Paso_S Line 2633  void Rectangle::assemblePDESystem(Paso_S
2633                          if (B.actsExpanded()) {                          if (B.actsExpanded()) {
2634                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2635                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2636                                      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)];
2637                                      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)];
2638                                      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)];
2639                                      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)];
2640                                      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)];
2641                                      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)];
2642                                      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)];
2643                                      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)];
2644                                      const register double tmp0_0 = B_1_0 + B_1_1;                                      const double tmp0_0 = B_1_0 + B_1_1;
2645                                      const register double tmp1_0 = B_1_2 + B_1_3;                                      const double tmp1_0 = B_1_2 + B_1_3;
2646                                      const register double tmp2_0 = B_0_1 + B_0_3;                                      const double tmp2_0 = B_0_1 + B_0_3;
2647                                      const register double tmp3_0 = B_0_0 + B_0_2;                                      const double tmp3_0 = B_0_0 + B_0_2;
2648                                      const register double tmp63_1 = B_1_1*w42;                                      const double tmp63_1 = B_1_1*w42;
2649                                      const register double tmp79_1 = B_1_1*w40;                                      const double tmp79_1 = B_1_1*w40;
2650                                      const register double tmp37_1 = tmp3_0*w35;                                      const double tmp37_1 = tmp3_0*w35;
2651                                      const register double tmp8_1 = tmp0_0*w32;                                      const double tmp8_1 = tmp0_0*w32;
2652                                      const register double tmp71_1 = B_0_1*w34;                                      const double tmp71_1 = B_0_1*w34;
2653                                      const register double tmp19_1 = B_0_3*w31;                                      const double tmp19_1 = B_0_3*w31;
2654                                      const register double tmp15_1 = B_0_3*w34;                                      const double tmp15_1 = B_0_3*w34;
2655                                      const register double tmp9_1 = tmp3_0*w34;                                      const double tmp9_1 = tmp3_0*w34;
2656                                      const register double tmp35_1 = B_1_0*w36;                                      const double tmp35_1 = B_1_0*w36;
2657                                      const register double tmp66_1 = B_0_3*w28;                                      const double tmp66_1 = B_0_3*w28;
2658                                      const register double tmp28_1 = B_1_0*w42;                                      const double tmp28_1 = B_1_0*w42;
2659                                      const register double tmp22_1 = B_1_0*w40;                                      const double tmp22_1 = B_1_0*w40;
2660                                      const register double tmp16_1 = B_1_2*w29;                                      const double tmp16_1 = B_1_2*w29;
2661                                      const register double tmp6_1 = tmp2_0*w35;                                      const double tmp6_1 = tmp2_0*w35;
2662                                      const register double tmp55_1 = B_1_3*w40;                                      const double tmp55_1 = B_1_3*w40;
2663                                      const register double tmp50_1 = B_1_3*w42;                                      const double tmp50_1 = B_1_3*w42;
2664                                      const register double tmp7_1 = tmp1_0*w29;                                      const double tmp7_1 = tmp1_0*w29;
2665                                      const register double tmp1_1 = tmp1_0*w32;                                      const double tmp1_1 = tmp1_0*w32;
2666                                      const register double tmp57_1 = B_0_3*w30;                                      const double tmp57_1 = B_0_3*w30;
2667                                      const register double tmp18_1 = B_1_1*w32;                                      const double tmp18_1 = B_1_1*w32;
2668                                      const register double tmp53_1 = B_1_0*w41;                                      const double tmp53_1 = B_1_0*w41;
2669                                      const register double tmp61_1 = B_1_3*w36;                                      const double tmp61_1 = B_1_3*w36;
2670                                      const register double tmp27_1 = B_0_3*w38;                                      const double tmp27_1 = B_0_3*w38;
2671                                      const register double tmp64_1 = B_0_2*w30;                                      const double tmp64_1 = B_0_2*w30;
2672                                      const register double tmp76_1 = B_0_1*w38;                                      const double tmp76_1 = B_0_1*w38;
2673                                      const register double tmp39_1 = tmp2_0*w34;                                      const double tmp39_1 = tmp2_0*w34;
2674                                      const register double tmp62_1 = B_0_1*w31;                                      const double tmp62_1 = B_0_1*w31;
2675                                      const register double tmp56_1 = B_0_0*w31;                                      const double tmp56_1 = B_0_0*w31;
2676                                      const register double tmp49_1 = B_1_1*w36;                                      const double tmp49_1 = B_1_1*w36;
2677                                      const register double tmp2_1 = B_0_2*w31;                                      const double tmp2_1 = B_0_2*w31;
2678                                      const register double tmp23_1 = B_0_2*w33;                                      const double tmp23_1 = B_0_2*w33;
2679                                      const register double tmp38_1 = B_1_1*w43;                                      const double tmp38_1 = B_1_1*w43;
2680                                      const register double tmp74_1 = B_1_2*w41;                                      const double tmp74_1 = B_1_2*w41;
2681                                      const register double tmp43_1 = B_1_1*w41;                                      const double tmp43_1 = B_1_1*w41;
2682                                      const register double tmp58_1 = B_0_2*w28;                                      const double tmp58_1 = B_0_2*w28;
2683                                      const register double tmp67_1 = B_0_0*w33;                                      const double tmp67_1 = B_0_0*w33;
2684                                      const register double tmp33_1 = tmp0_0*w39;                                      const double tmp33_1 = tmp0_0*w39;
2685                                      const register double tmp4_1 = B_0_0*w28;                                      const double tmp4_1 = B_0_0*w28;
2686                                      const register double tmp20_1 = B_0_0*w30;                                      const double tmp20_1 = B_0_0*w30;
2687                                      const register double tmp13_1 = B_0_2*w38;                                      const double tmp13_1 = B_0_2*w38;
2688                                      const register double tmp65_1 = B_1_2*w43;                                      const double tmp65_1 = B_1_2*w43;
2689                                      const register double tmp0_1 = tmp0_0*w29;                                      const double tmp0_1 = tmp0_0*w29;
2690                                      const register double tmp41_1 = tmp3_0*w33;                                      const double tmp41_1 = tmp3_0*w33;
2691                                      const register double tmp73_1 = B_0_2*w37;                                      const double tmp73_1 = B_0_2*w37;
2692                                      const register double tmp69_1 = B_0_0*w38;                                      const double tmp69_1 = B_0_0*w38;
2693                                      const register double tmp48_1 = B_1_2*w39;                                      const double tmp48_1 = B_1_2*w39;
2694                                      const register double tmp59_1 = B_0_1*w33;                                      const double tmp59_1 = B_0_1*w33;
2695                                      const register double tmp17_1 = B_1_3*w41;                                      const double tmp17_1 = B_1_3*w41;
2696                                      const register double tmp5_1 = B_0_3*w33;                                      const double tmp5_1 = B_0_3*w33;
2697                                      const register double tmp3_1 = B_0_1*w30;                                      const double tmp3_1 = B_0_1*w30;
2698                                      const register double tmp21_1 = B_0_1*w28;                                      const double tmp21_1 = B_0_1*w28;
2699                                      const register double tmp42_1 = B_1_0*w29;                                      const double tmp42_1 = B_1_0*w29;
2700                                      const register double tmp54_1 = B_1_2*w32;                                      const double tmp54_1 = B_1_2*w32;
2701                                      const register double tmp60_1 = B_1_0*w39;                                      const double tmp60_1 = B_1_0*w39;
2702                                      const register double tmp32_1 = tmp1_0*w36;                                      const double tmp32_1 = tmp1_0*w36;
2703                                      const register double tmp10_1 = B_0_1*w37;                                      const double tmp10_1 = B_0_1*w37;
2704                                      const register double tmp14_1 = B_0_0*w35;                                      const double tmp14_1 = B_0_0*w35;
2705                                      const register double tmp29_1 = B_0_1*w35;                                      const double tmp29_1 = B_0_1*w35;
2706                                      const register double tmp26_1 = B_1_2*w36;                                      const double tmp26_1 = B_1_2*w36;
2707                                      const register double tmp30_1 = B_1_3*w43;                                      const double tmp30_1 = B_1_3*w43;
2708                                      const register double tmp70_1 = B_0_2*w35;                                      const double tmp70_1 = B_0_2*w35;
2709                                      const register double tmp34_1 = B_1_3*w39;                                      const double tmp34_1 = B_1_3*w39;
2710                                      const register double tmp51_1 = B_1_0*w43;                                      const double tmp51_1 = B_1_0*w43;
2711                                      const register double tmp31_1 = B_0_2*w34;                                      const double tmp31_1 = B_0_2*w34;
2712                                      const register double tmp45_1 = tmp3_0*w28;                                      const double tmp45_1 = tmp3_0*w28;
2713                                      const register double tmp11_1 = tmp1_0*w39;                                      const double tmp11_1 = tmp1_0*w39;
2714                                      const register double tmp52_1 = B_1_1*w29;                                      const double tmp52_1 = B_1_1*w29;
2715                                      const register double tmp44_1 = B_1_3*w32;                                      const double tmp44_1 = B_1_3*w32;
2716                                      const register double tmp25_1 = B_1_1*w39;                                      const double tmp25_1 = B_1_1*w39;
2717                                      const register double tmp47_1 = tmp2_0*w33;                                      const double tmp47_1 = tmp2_0*w33;
2718                                      const register double tmp72_1 = B_1_3*w29;                                      const double tmp72_1 = B_1_3*w29;
2719                                      const register double tmp40_1 = tmp2_0*w28;                                      const double tmp40_1 = tmp2_0*w28;
2720                                      const register double tmp46_1 = B_1_2*w40;                                      const double tmp46_1 = B_1_2*w40;
2721                                      const register double tmp36_1 = B_1_2*w42;                                      const double tmp36_1 = B_1_2*w42;
2722                                      const register double tmp24_1 = B_0_0*w37;                                      const double tmp24_1 = B_0_0*w37;
2723                                      const register double tmp77_1 = B_0_3*w35;                                      const double tmp77_1 = B_0_3*w35;
2724                                      const register double tmp68_1 = B_0_3*w37;                                      const double tmp68_1 = B_0_3*w37;
2725                                      const register double tmp78_1 = B_0_0*w34;                                      const double tmp78_1 = B_0_0*w34;
2726                                      const register double tmp12_1 = tmp0_0*w36;                                      const double tmp12_1 = tmp0_0*w36;
2727                                      const register double tmp75_1 = B_1_0*w32;                                      const double tmp75_1 = B_1_0*w32;
2728                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2729                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2730                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
# Line 2768  void Rectangle::assemblePDESystem(Paso_S Line 2746  void Rectangle::assemblePDESystem(Paso_S
2746                          } else { // constant data                          } else { // constant data
2747                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2748                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2749                                      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)];
2750                                      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)];
2751                                      const register double tmp6_1 = B_1*w50;                                      const double tmp6_1 = B_1*w50;
2752                                      const register double tmp1_1 = B_1*w45;                                      const double tmp1_1 = B_1*w45;
2753                                      const register double tmp5_1 = B_1*w49;                                      const double tmp5_1 = B_1*w49;
2754                                      const register double tmp4_1 = B_1*w48;                                      const double tmp4_1 = B_1*w48;
2755                                      const register double tmp0_1 = B_0*w44;                                      const double tmp0_1 = B_0*w44;
2756                                      const register double tmp2_1 = B_0*w46;                                      const double tmp2_1 = B_0*w46;
2757                                      const register double tmp7_1 = B_0*w51;                                      const double tmp7_1 = B_0*w51;
2758                                      const register double tmp3_1 = B_0*w47;                                      const double tmp3_1 = B_0*w47;
2759                                      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;
2760                                      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;
2761                                      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 2807  void Rectangle::assemblePDESystem(Paso_S Line 2785  void Rectangle::assemblePDESystem(Paso_S
2785                          if (C.actsExpanded()) {                          if (C.actsExpanded()) {
2786                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2787                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2788                                      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)];
2789                                      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)];
2790                                      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)];
2791                                      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)];
2792                                      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)];
2793                                      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)];
2794                                      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)];
2795                                      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)];
2796                                      const register double tmp2_0 = C_0_1 + C_0_3;                                      const double tmp2_0 = C_0_1 + C_0_3;
2797                                      const register double tmp1_0 = C_1_2 + C_1_3;                                      const double tmp1_0 = C_1_2 + C_1_3;
2798                                      const register double tmp3_0 = C_0_0 + C_0_2;                                      const double tmp3_0 = C_0_0 + C_0_2;
2799                                      const register double tmp0_0 = C_1_0 + C_1_1;                                      const double tmp0_0 = C_1_0 + C_1_1;
2800                                      const register double tmp64_1 = C_0_2*w30;                                      const double tmp64_1 = C_0_2*w30;
2801                                      const register double tmp14_1 = C_0_2*w28;                                      const double tmp14_1 = C_0_2*w28;
2802                                      const register double tmp19_1 = C_0_3*w31;                                      const double tmp19_1 = C_0_3*w31;
2803                                      const register double tmp22_1 = C_1_0*w40;                                      const double tmp22_1 = C_1_0*w40;
2804                                      const register double tmp37_1 = tmp3_0*w35;                                      const double tmp37_1 = tmp3_0*w35;
2805                                      const register double tmp29_1 = C_0_1*w35;                                      const double tmp29_1 = C_0_1*w35;
2806                                      const register double tmp73_1 = C_0_2*w37;                                      const double tmp73_1 = C_0_2*w37;
2807                                      const register double tmp74_1 = C_1_2*w41;                                      const double tmp74_1 = C_1_2*w41;
2808                                      const register double tmp52_1 = C_1_3*w39;                                      const double tmp52_1 = C_1_3*w39;
2809                                      const register double tmp25_1 = C_1_1*w39;                                      const double tmp25_1 = C_1_1*w39;
2810                                      const register double tmp62_1 = C_0_1*w31;                                      const double tmp62_1 = C_0_1*w31;
2811                                      const register double tmp79_1 = C_1_1*w40;                                      const double tmp79_1 = C_1_1*w40;
2812                                      const register double tmp43_1 = C_1_1*w36;                                      const double tmp43_1 = C_1_1*w36;
2813                                      const register double tmp27_1 = C_0_3*w38;                                      const double tmp27_1 = C_0_3*w38;
2814                                      const register double tmp28_1 = C_1_0*w42;                                      const double tmp28_1 = C_1_0*w42;
2815                                      const register double tmp63_1 = C_1_1*w42;                                      const double tmp63_1 = C_1_1*w42;
2816                                      const register double tmp59_1 = C_0_3*w34;                                      const double tmp59_1 = C_0_3*w34;
2817                                      const register double tmp72_1 = C_1_3*w29;                                      const double tmp72_1 = C_1_3*w29;
2818                                      const register double tmp40_1 = tmp2_0*w35;                                      const double tmp40_1 = tmp2_0*w35;
2819                                      const register double tmp13_1 = C_0_3*w30;                                      const double tmp13_1 = C_0_3*w30;
2820                                      const register double tmp51_1 = C_1_2*w40;                                      const double tmp51_1 = C_1_2*w40;
2821                                      const register double tmp54_1 = C_1_2*w42;                                      const double tmp54_1 = C_1_2*w42;
2822                                      const register double tmp12_1 = C_0_0*w31;                                      const double tmp12_1 = C_0_0*w31;
2823                                      const register double tmp2_1 = tmp1_0*w32;                                      const double tmp2_1 = tmp1_0*w32;
2824                                      const register double tmp68_1 = C_0_2*w31;                                      const double tmp68_1 = C_0_2*w31;
2825                                      const register double tmp75_1 = C_1_0*w32;                                      const double tmp75_1 = C_1_0*w32;
2826                                      const register double tmp49_1 = C_1_1*w41;                                      const double tmp49_1 = C_1_1*w41;
2827                                      const register double tmp4_1 = C_0_2*w35;                                      const double tmp4_1 = C_0_2*w35;
2828                                      const register double tmp66_1 = C_0_3*w28;                                      const double tmp66_1 = C_0_3*w28;
2829                                      const register double tmp56_1 = C_0_1*w37;                                      const double tmp56_1 = C_0_1*w37;
2830                                      const register double tmp5_1 = C_0_1*w34;                                      const double tmp5_1 = C_0_1*w34;
2831                                      const register double tmp38_1 = tmp2_0*w34;                                      const double tmp38_1 = tmp2_0*w34;
2832                                      const register double tmp76_1 = C_0_1*w38;                                      const double tmp76_1 = C_0_1*w38;
2833                                      const register double tmp21_1 = C_0_1*w28;                                      const double tmp21_1 = C_0_1*w28;
2834                                      const register double tmp69_1 = C_0_1*w30;                                      const double tmp69_1 = C_0_1*w30;
2835                                      const register double tmp53_1 = C_1_0*w36;                                      const double tmp53_1 = C_1_0*w36;
2836                                      const register double tmp42_1 = C_1_2*w39;                                      const double tmp42_1 = C_1_2*w39;
2837                                      const register double tmp32_1 = tmp1_0*w29;                                      const double tmp32_1 = tmp1_0*w29;
2838                                      const register double tmp45_1 = C_1_0*w43;                                      const double tmp45_1 = C_1_0*w43;
2839                                      const register double tmp33_1 = tmp0_0*w32;                                      const double tmp33_1 = tmp0_0*w32;
2840                                      const register double tmp35_1 = C_1_0*w41;                                      const double tmp35_1 = C_1_0*w41;
2841                                      const register double tmp26_1 = C_1_2*w36;                                      const double tmp26_1 = C_1_2*w36;
2842                                      const register double tmp67_1 = C_0_0*w33;                                      const double tmp67_1 = C_0_0*w33;
2843                                      const register double tmp31_1 = C_0_2*w34;                                      const double tmp31_1 = C_0_2*w34;
2844                                      const register double tmp20_1 = C_0_0*w30;                                      const double tmp20_1 = C_0_0*w30;
2845                                      const register double tmp70_1 = C_0_0*w28;                                      const double tmp70_1 = C_0_0*w28;
2846                                      const register double tmp8_1 = tmp0_0*w39;                                      const double tmp8_1 = tmp0_0*w39;
2847                                      const register double tmp30_1 = C_1_3*w43;                                      const double tmp30_1 = C_1_3*w43;
2848                                      const register double tmp0_1 = tmp0_0*w29;                                      const double tmp0_1 = tmp0_0*w29;
2849                                      const register double tmp17_1 = C_1_3*w41;                                      const double tmp17_1 = C_1_3*w41;
2850                                      const register double tmp58_1 = C_0_0*w35;                                      const double tmp58_1 = C_0_0*w35;
2851                                      const register double tmp9_1 = tmp3_0*w33;                                      const double tmp9_1 = tmp3_0*w33;
2852                                      const register double tmp61_1 = C_1_3*w36;                                      const double tmp61_1 = C_1_3*w36;
2853                                      const register double tmp41_1 = tmp3_0*w34;                                      const double tmp41_1 = tmp3_0*w34;
2854                                      const register double tmp50_1 = C_1_3*w32;                                      const double tmp50_1 = C_1_3*w32;
2855                                      const register double tmp18_1 = C_1_1*w32;                                      const double tmp18_1 = C_1_1*w32;
2856                                      const register double tmp6_1 = tmp1_0*w36;                                      const double tmp6_1 = tmp1_0*w36;
2857                                      const register double tmp3_1 = C_0_0*w38;                                      const double tmp3_1 = C_0_0*w38;
2858                                      const register double tmp34_1 = C_1_1*w29;                                      const double tmp34_1 = C_1_1*w29;
2859                                      const register double tmp77_1 = C_0_3*w35;                                      const double tmp77_1 = C_0_3*w35;
2860                                      const register double tmp65_1 = C_1_2*w43;                                      const double tmp65_1 = C_1_2*w43;
2861                                      const register double tmp71_1 = C_0_3*w33;                                      const double tmp71_1 = C_0_3*w33;
2862                                      const register double tmp55_1 = C_1_1*w43;                                      const double tmp55_1 = C_1_1*w43;
2863                                      const register double tmp46_1 = tmp3_0*w28;                                      const double tmp46_1 = tmp3_0*w28;
2864                                      const register double tmp24_1 = C_0_0*w37;                                      const double tmp24_1 = C_0_0*w37;
2865                                      const register double tmp10_1 = tmp1_0*w39;                                      const double tmp10_1 = tmp1_0*w39;
2866                                      const register double tmp48_1 = C_1_0*w29;                                      const double tmp48_1 = C_1_0*w29;
2867                                      const register double tmp15_1 = C_0_1*w33;                                      const double tmp15_1 = C_0_1*w33;
2868                                      const register double tmp36_1 = C_1_2*w32;                                      const double tmp36_1 = C_1_2*w32;
2869                                      const register double tmp60_1 = C_1_0*w39;                                      const double tmp60_1 = C_1_0*w39;
2870                                      const register double tmp47_1 = tmp2_0*w33;                                      const double tmp47_1 = tmp2_0*w33;
2871                                      const register double tmp16_1 = C_1_2*w29;                                      const double tmp16_1 = C_1_2*w29;
2872                                      const register double tmp1_1 = C_0_3*w37;                                      const double tmp1_1 = C_0_3*w37;
2873                                      const register double tmp7_1 = tmp2_0*w28;                                      const double tmp7_1 = tmp2_0*w28;
2874                                      const register double tmp39_1 = C_1_3*w40;                                      const double tmp39_1 = C_1_3*w40;
2875                                      const register double tmp44_1 = C_1_3*w42;                                      const double tmp44_1 = C_1_3*w42;
2876                                      const register double tmp57_1 = C_0_2*w38;                                      const double tmp57_1 = C_0_2*w38;
2877                                      const register double tmp78_1 = C_0_0*w34;                                      const double tmp78_1 = C_0_0*w34;
2878                                      const register double tmp11_1 = tmp0_0*w36;                                      const double tmp11_1 = tmp0_0*w36;
2879                                      const register double tmp23_1 = C_0_2*w33;                                      const double tmp23_1 = C_0_2*w33;
2880                                      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;
2881                                      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;
2882                                      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 2920  void Rectangle::assemblePDESystem(Paso_S Line 2898  void Rectangle::assemblePDESystem(Paso_S
2898                          } else { // constant data                          } else { // constant data
2899                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2900                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2901                                      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)];
2902                                      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)];
2903                                      const register double tmp1_1 = C_1*w45;                                      const double tmp1_1 = C_1*w45;
2904                                      const register double tmp3_1 = C_0*w51;                                      const double tmp3_1 = C_0*w51;
2905                                      const register double tmp4_1 = C_0*w44;                                      const double tmp4_1 = C_0*w44;
2906                                      const register double tmp7_1 = C_0*w46;                                      const double tmp7_1 = C_0*w46;
2907                                      const register double tmp5_1 = C_1*w49;                                      const double tmp5_1 = C_1*w49;
2908                                      const register double tmp2_1 = C_1*w48;                                      const double tmp2_1 = C_1*w48;
2909                                      const register double tmp0_1 = C_0*w47;                                      const double tmp0_1 = C_0*w47;
2910                                      const register double tmp6_1 = C_1*w50;                                      const double tmp6_1 = C_1*w50;
2911                                      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;
2912                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
2913                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
# Line 2959  void Rectangle::assemblePDESystem(Paso_S Line 2937  void Rectangle::assemblePDESystem(Paso_S
2937                          if (D.actsExpanded()) {                          if (D.actsExpanded()) {
2938                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2939                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2940                                      const register double D_0 = D_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double D_0 = D_p[INDEX3(k, m, 0, numEq, numComp)];
2941                                      const register double D_1 = D_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double D_1 = D_p[INDEX3(k, m, 1, numEq, numComp)];
2942                                      const register double D_2 = D_p[INDEX3(k, m, 2, numEq, numComp)];                                      const double D_2 = D_p[INDEX3(k, m, 2, numEq, numComp)];
2943                                      const register double D_3 = D_p[INDEX3(k, m, 3, numEq, numComp)];                                      const double D_3 = D_p[INDEX3(k, m, 3, numEq, numComp)];
2944                                      const register double tmp4_0 = D_1 + D_3;                                      const double tmp4_0 = D_1 + D_3;
2945                                      const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;                                      const double tmp2_0 = D_0 + D_1 + D_2 + D_3;
2946                                      const register double tmp5_0 = D_0 + D_2;                                      const double tmp5_0 = D_0 + D_2;
2947                                      const register double tmp0_0 = D_0 + D_1;                                      const double tmp0_0 = D_0 + D_1;
2948                                      const register double tmp6_0 = D_0 + D_3;                                      const double tmp6_0 = D_0 + D_3;
2949                                      const register double tmp1_0 = D_2 + D_3;                                      const double tmp1_0 = D_2 + D_3;
2950                                      const register double tmp3_0 = D_1 + D_2;                                      const double tmp3_0 = D_1 + D_2;
2951                                      const register double tmp16_1 = D_1*w56;                                      const double tmp16_1 = D_1*w56;
2952                                      const register double tmp14_1 = tmp6_0*w54;                                      const double tmp14_1 = tmp6_0*w54;
2953                                      const register double tmp8_1 = D_3*w55;                                      const double tmp8_1 = D_3*w55;
2954                                      const register double tmp2_1 = tmp2_0*w54;                                      const double tmp2_1 = tmp2_0*w54;
2955                                      const register double tmp12_1 = tmp5_0*w52;                                      const double tmp12_1 = tmp5_0*w52;
2956                                      const register double tmp4_1 = tmp0_0*w53;                                      const double tmp4_1 = tmp0_0*w53;
2957                                      const register double tmp3_1 = tmp1_0*w52;                                      const double tmp3_1 = tmp1_0*w52;
2958                                      const register double tmp13_1 = tmp4_0*w53;                                      const double tmp13_1 = tmp4_0*w53;
2959                                      const register double tmp10_1 = tmp4_0*w52;                                      const double tmp10_1 = tmp4_0*w52;
2960                                      const register double tmp1_1 = tmp1_0*w53;                                      const double tmp1_1 = tmp1_0*w53;
2961                                      const register double tmp7_1 = D_3*w56;                                      const double tmp7_1 = D_3*w56;
2962                                      const register double tmp0_1 = tmp0_0*w52;                                      const double tmp0_1 = tmp0_0*w52;
2963                                      const register double tmp11_1 = tmp5_0*w53;                                      const double tmp11_1 = tmp5_0*w53;
2964                                      const register double tmp9_1 = D_0*w56;                                      const double tmp9_1 = D_0*w56;
2965                                      const register double tmp5_1 = tmp3_0*w54;                                      const double tmp5_1 = tmp3_0*w54;
2966                                      const register double tmp18_1 = D_2*w56;                                      const double tmp18_1 = D_2*w56;
2967                                      const register double tmp17_1 = D_1*w55;                                      const double tmp17_1 = D_1*w55;
2968                                      const register double tmp6_1 = D_0*w55;                                      const double tmp6_1 = D_0*w55;
2969                                      const register double tmp15_1 = D_2*w55;                                      const double tmp15_1 = D_2*w55;
2970                                      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;
2971                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1;
2972                                      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 3010  void Rectangle::assemblePDESystem(Paso_S Line 2988  void Rectangle::assemblePDESystem(Paso_S
2988                          } else { // constant data                          } else { // constant data
2989                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2990                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2991                                      const register double D_0 = D_p[INDEX2(k, m, numEq)];                                      const double D_0 = D_p[INDEX2(k, m, numEq)];
2992                                      const register double tmp2_1 = D_0*w59;                                      const double tmp2_1 = D_0*w59;
2993                                      const register double tmp1_1 = D_0*w58;                                      const double tmp1_1 = D_0*w58;
2994                                      const register double tmp0_1 = D_0*w57;                                      const double tmp0_1 = D_0*w57;
2995                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
2996                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1;
2997                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
# Line 3042  void Rectangle::assemblePDESystem(Paso_S Line 3020  void Rectangle::assemblePDESystem(Paso_S
3020                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3021                          if (X.actsExpanded()) {                          if (X.actsExpanded()) {
3022                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3023                                  const register double X_0_0 = X_p[INDEX3(k, 0, 0, numEq, 2)];                                  const double X_0_0 = X_p[INDEX3(k, 0, 0, numEq, 2)];
3024                                  const register double X_1_0 = X_p[INDEX3(k, 1, 0, numEq, 2)];                                  const double X_1_0 = X_p[INDEX3(k, 1, 0, numEq, 2)];
3025                                  const register double X_0_1 = X_p[INDEX3(k, 0, 1, numEq, 2)];                                  const double X_0_1 = X_p[INDEX3(k, 0, 1, numEq, 2)];
3026                                  const register double X_1_1 = X_p[INDEX3(k, 1, 1, numEq, 2)];                                  const double X_1_1 = X_p[INDEX3(k, 1, 1, numEq, 2)];
3027                                  const register double X_0_2 = X_p[INDEX3(k, 0, 2, numEq, 2)];                                  const double X_0_2 = X_p[INDEX3(k, 0, 2, numEq, 2)];
3028                                  const register double X_1_2 = X_p[INDEX3(k, 1, 2, numEq, 2)];                                  const double X_1_2 = X_p[INDEX3(k, 1, 2, numEq, 2)];
3029                                  const register double X_0_3 = X_p[INDEX3(k, 0, 3, numEq, 2)];                                  const double X_0_3 = X_p[INDEX3(k, 0, 3, numEq, 2)];
3030                                  const register double X_1_3 = X_p[INDEX3(k, 1, 3, numEq, 2)];                                  const double X_1_3 = X_p[INDEX3(k, 1, 3, numEq, 2)];
3031                                  const register double tmp1_0 = X_1_1 + X_1_3;                                  const double tmp1_0 = X_1_1 + X_1_3;
3032                                  const register double tmp3_0 = X_0_0 + X_0_1;                                  const double tmp3_0 = X_0_0 + X_0_1;
3033                                  const register double tmp2_0 = X_1_0 + X_1_2;                                  const double tmp2_0 = X_1_0 + X_1_2;
3034                                  const register double tmp0_0 = X_0_2 + X_0_3;                                  const double tmp0_0 = X_0_2 + X_0_3;
3035                                  const register double tmp8_1 = tmp2_0*w66;                                  const double tmp8_1 = tmp2_0*w66;
3036                                  const register double tmp5_1 = tmp3_0*w64;                                  const double tmp5_1 = tmp3_0*w64;
3037                                  const register double tmp14_1 = tmp0_0*w64;                                  const double tmp14_1 = tmp0_0*w64;
3038                                  const register double tmp3_1 = tmp3_0*w60;                                  const double tmp3_1 = tmp3_0*w60;
3039                                  const register double tmp9_1 = tmp3_0*w63;                                  const double tmp9_1 = tmp3_0*w63;
3040                                  const register double tmp13_1 = tmp3_0*w65;                                  const double tmp13_1 = tmp3_0*w65;
3041                                  const register double tmp12_1 = tmp1_0*w66;                                  const double tmp12_1 = tmp1_0*w66;
3042                                  const register double tmp10_1 = tmp0_0*w60;                                  const double tmp10_1 = tmp0_0*w60;
3043                                  const register double tmp2_1 = tmp2_0*w61;                                  const double tmp2_1 = tmp2_0*w61;
3044                                  const register double tmp6_1 = tmp2_0*w62;                                  const double tmp6_1 = tmp2_0*w62;
3045                                  const register double tmp4_1 = tmp0_0*w65;                                  const double tmp4_1 = tmp0_0*w65;
3046                                  const register double tmp11_1 = tmp1_0*w67;                                  const double tmp11_1 = tmp1_0*w67;
3047                                  const register double tmp1_1 = tmp1_0*w62;                                  const double tmp1_1 = tmp1_0*w62;
3048                                  const register double tmp7_1 = tmp1_0*w61;                                  const double tmp7_1 = tmp1_0*w61;
3049                                  const register double tmp0_1 = tmp0_0*w63;                                  const double tmp0_1 = tmp0_0*w63;
3050                                  const register double tmp15_1 = tmp2_0*w67;                                  const double tmp15_1 = tmp2_0*w67;
3051                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3052                                  EM_F[INDEX2(k,1,numEq)]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
3053                                  EM_F[INDEX2(k,2,numEq)]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
# Line 3077  void Rectangle::assemblePDESystem(Paso_S Line 3055  void Rectangle::assemblePDESystem(Paso_S
3055                              }                              }
3056                          } else { // constant data                          } else { // constant data
3057                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3058                                  const register double X_0 = X_p[INDEX2(k, 0, numEq)];                                  const double X_0 = X_p[INDEX2(k, 0, numEq)];
3059                                  const register double X_1 = X_p[INDEX2(k, 1, numEq)];                                  const double X_1 = X_p[INDEX2(k, 1, numEq)];
3060                                  const register double tmp0_1 = X_1*w69;                                  const double tmp0_1 = X_1*w69;
3061                                  const register double tmp1_1 = X_0*w68;                                  const double tmp1_1 = X_0*w68;
3062                                  const register double tmp2_1 = X_0*w70;                                  const double tmp2_1 = X_0*w70;
3063                                  const register double tmp3_1 = X_1*w71;                                  const double tmp3_1 = X_1*w71;
3064                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3065                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;
3066                                  EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;
# Line 3098  void Rectangle::assemblePDESystem(Paso_S Line 3076  void Rectangle::assemblePDESystem(Paso_S
3076                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3077                          if (Y.actsExpanded()) {                          if (Y.actsExpanded()) {
3078                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3079                                  const register double Y_0 = Y_p[INDEX2(k, 0, numEq)];                                  const double Y_0 = Y_p[INDEX2(k, 0, numEq)];
3080                                  const register double Y_1 = Y_p[INDEX2(k, 1, numEq)];                                  const double Y_1 = Y_p[INDEX2(k, 1, numEq)];
3081                                  const register double Y_2 = Y_p[INDEX2(k, 2, numEq)];                                  const double Y_2 = Y_p[INDEX2(k, 2, numEq)];
3082                                  const register double Y_3 = Y_p[INDEX2(k, 3, numEq)];                                  const double Y_3 = Y_p[INDEX2(k, 3, numEq)];
3083                                  const register double tmp0_0 = Y_1 + Y_2;                                  const double tmp0_0 = Y_1 + Y_2;
3084                                  const register double tmp1_0 = Y_0 + Y_3;                                  const double tmp1_0 = Y_0 + Y_3;
3085                                  const register double tmp0_1 = Y_0*w72;                                  const double tmp0_1 = Y_0*w72;
3086                                  const register double tmp1_1 = tmp0_0*w73;                                  const double tmp1_1 = tmp0_0*w73;
3087                                  const register double tmp2_1 = Y_3*w74;                                  const double tmp2_1 = Y_3*w74;
3088                                  const register double tmp3_1 = Y_1*w72;                                  const double tmp3_1 = Y_1*w72;
3089                                  const register double tmp4_1 = tmp1_0*w73;                                  const double tmp4_1 = tmp1_0*w73;
3090                                  const register double tmp5_1 = Y_2*w74;                                  const double tmp5_1 = Y_2*w74;
3091                                  const register double tmp6_1 = Y_2*w72;                                  const double tmp6_1 = Y_2*w72;
3092                                  const register double tmp7_1 = Y_1*w74;                                  const double tmp7_1 = Y_1*w74;
3093                                  const register double tmp8_1 = Y_3*w72;                                  const double tmp8_1 = Y_3*w72;
3094                                  const register double tmp9_1 = Y_0*w74;                                  const double tmp9_1 = Y_0*w74;
3095                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
3096                                  EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
3097                                  EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
# Line 3121  void Rectangle::assemblePDESystem(Paso_S Line 3099  void Rectangle::assemblePDESystem(Paso_S
3099                              }                              }
3100                          } else { // constant data                          } else { // constant data
3101                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3102                                  const register double tmp0_1 = Y_p[k]*w75;                                  const double tmp0_1 = Y_p[k]*w75;
3103                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3104                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3105                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
# Line 3132  void Rectangle::assemblePDESystem(Paso_S Line 3110  void Rectangle::assemblePDESystem(Paso_S
3110    
3111                      // 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)
3112                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
3113                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
3114                      rowIndex.push_back(m_dofMap[firstNode]);                              firstNode, numEq, numComp);
                     rowIndex.push_back(m_dofMap[firstNode+1]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0]);  
                     rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);  
                     if (add_EM_F) {  
                         double *F_p=rhs.getSampleDataRW(0);  
                         for (index_t i=0; i<rowIndex.size(); i++) {  
                             if (rowIndex[i]<getNumDOF()) {  
                                 for (index_t eq=0; eq<numEq; eq++) {  
                                     F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];  
                                 }  
                             }  
                         }  
                     }  
                     if (add_EM_S) {  
                         addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);  
                     }  
   
3115                  } // end k0 loop                  } // end k0 loop
3116              } // end k1 loop              } // end k1 loop
3117          } // end of colouring          } // end of colouring
# Line 3161  void Rectangle::assemblePDESystem(Paso_S Line 3122  void Rectangle::assemblePDESystem(Paso_S
3122  void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,
3123          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
3124          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
3125          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  
3126  {  {
3127      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
3128      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 3211  void Rectangle::assemblePDESystemReduced Line 3171  void Rectangle::assemblePDESystemReduced
3171                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
3172                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3173                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3174                                  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)];
3175                                  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)];
3176                                  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)];
3177                                  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)];
3178                                  const register double tmp0_0 = A_01 + A_10;                                  const double tmp0_0 = A_01 + A_10;
3179                                  const register double tmp0_1 = A_11*w3;                                  const double tmp0_1 = A_11*w3;
3180                                  const register double tmp1_1 = A_00*w0;                                  const double tmp1_1 = A_00*w0;
3181                                  const register double tmp2_1 = A_01*w1;                                  const double tmp2_1 = A_01*w1;
3182                                  const register double tmp3_1 = A_10*w2;                                  const double tmp3_1 = A_10*w2;
3183                                  const register double tmp4_1 = tmp0_0*w1;                                  const double tmp4_1 = tmp0_0*w1;
3184                                  const register double tmp5_1 = A_11*w4;                                  const double tmp5_1 = A_11*w4;
3185                                  const register double tmp6_1 = A_00*w5;                                  const double tmp6_1 = A_00*w5;
3186                                  const register double tmp7_1 = tmp0_0*w2;                                  const double tmp7_1 = tmp0_0*w2;
3187                                  const register double tmp8_1 = A_10*w1;                                  const double tmp8_1 = A_10*w1;
3188                                  const register double tmp9_1 = A_01*w2;                                  const double tmp9_1 = A_01*w2;
3189                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
3190                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
3191                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
# Line 3253  void Rectangle::assemblePDESystemReduced Line 3213  void Rectangle::assemblePDESystemReduced
3213                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
3214                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3215                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3216                                  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)];
3217                                  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)];
3218                                  const register double tmp0_1 = B_0*w6;                                  const double tmp0_1 = B_0*w6;
3219                                  const register double tmp1_1 = B_1*w7;                                  const double tmp1_1 = B_1*w7;
3220                                  const register double tmp2_1 = B_0*w8;                                  const double tmp2_1 = B_0*w8;
3221                                  const register double tmp3_1 = B_1*w9;                                  const double tmp3_1 = B_1*w9;
3222                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3223                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3224                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
# Line 3286  void Rectangle::assemblePDESystemReduced Line 3246  void Rectangle::assemblePDESystemReduced
3246                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
3247                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3248                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3249                                  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)];
3250                                  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)];
3251                                  const register double tmp0_1 = C_0*w8;                                  const double tmp0_1 = C_0*w8;
3252                                  const register double tmp1_1 = C_1*w7;                                  const double tmp1_1 = C_1*w7;
3253                                  const register double tmp2_1 = C_1*w9;                                  const double tmp2_1 = C_1*w9;
3254                                  const register double tmp3_1 = C_0*w6;                                  const double tmp3_1 = C_0*w6;
3255                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3256                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3257                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
# Line 3319  void Rectangle::assemblePDESystemReduced Line 3279  void Rectangle::assemblePDESystemReduced
3279                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
3280                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3281                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3282                                  const register double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;                                  const double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;
3283                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
3284                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
3285                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;
# Line 3346  void Rectangle::assemblePDESystemReduced Line 3306  void Rectangle::assemblePDESystemReduced
3306                          add_EM_F=true;                          add_EM_F=true;
3307                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3308                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3309                              const register double X_0 = X_p[INDEX2(k, 0, numEq)];                              const double X_0 = X_p[INDEX2(k, 0, numEq)];
3310                              const register double X_1 = X_p[INDEX2(k, 1, numEq)];                              const double X_1 = X_p[INDEX2(k, 1, numEq)];
3311                              const register double tmp0_1 = X_0*w11;                              const double tmp0_1 = X_0*w11;
3312                              const register double tmp1_1 = X_1*w12;                              const double tmp1_1 = X_1*w12;
3313                              const register double tmp2_1 = X_0*w13;                              const double tmp2_1 = X_0*w13;
3314                              const register double tmp3_1 = X_1*w14;                              const double tmp3_1 = X_1*w14;
3315                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3316                              EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;
3317                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp3_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp3_1;
# Line 3365  void Rectangle::assemblePDESystemReduced Line 3325  void Rectangle::assemblePDESystemReduced
3325                          add_EM_F=true;                          add_EM_F=true;
3326                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3327                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3328                              const register double tmp0_1 = Y_p[k]*w15;                              const double tmp0_1 = Y_p[k]*w15;
3329                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3330                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3331                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
# Line 3375  void Rectangle::assemblePDESystemReduced Line 3335  void Rectangle::assemblePDESystemReduced
3335    
3336                      // 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)
3337                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
3338                      IndexVector rowIndex;                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F,
3339                      rowIndex.push_back(m_dofMap[firstNode]);                              firstNode, numEq, numComp);
3340                      rowIndex.push_back(m_dofMap[firstNode+1]);                  } // end k0 loop
3341                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);              } // end k1 loop
3342                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);          } // end of colouring
3343        } // end of parallel region
3344    }
3345    
3346    //protected
3347    void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
3348          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3349    {
3350        const double h0 = m_l0/m_gNE0;
3351        const double h1 = m_l1/m_gNE1;
3352        /* GENERATOR SNIP_PDEBC_SINGLE_PRE TOP */
3353        const double w0 = 0.31100423396407310779*h1;
3354        const double w1 = 0.022329099369260225539*h1;
3355        const double w10 = 0.022329099369260225539*h0;
3356        const double w11 = 0.16666666666666666667*h0;
3357        const double w12 = 0.33333333333333333333*h0;
3358        const double w13 = 0.39433756729740644113*h0;
3359        const double w14 = 0.10566243270259355887*h0;
3360        const double w15 = 0.5*h0;
3361        const double w2 = 0.083333333333333333333*h1;
3362        const double w3 = 0.33333333333333333333*h1;
3363        const double w4 = 0.16666666666666666667*h1;
3364        const double w5 = 0.39433756729740644113*h1;
3365        const double w6 = 0.10566243270259355887*h1;
3366        const double w7 = 0.5*h1;
3367        const double w8 = 0.083333333333333333333*h0;
3368        const double w9 = 0.31100423396407310779*h0;
3369        /* GENERATOR SNIP_PDEBC_SINGLE_PRE BOTTOM */
3370        const bool add_EM_S=!d.isEmpty();
3371        const bool add_EM_F=!y.isEmpty();
3372    #pragma omp parallel
3373        {
3374            if (m_faceOffset[0] > -1) {
3375                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3376    #pragma omp for nowait
3377                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3378                        vector<double> EM_S(4*4, 0);
3379                        vector<double> EM_F(4, 0);
3380                        const index_t e = k1;
3381                        /* GENERATOR SNIP_PDEBC_SINGLE_0 TOP */
3382                        ///////////////
3383                        // process d //
3384                        ///////////////
3385                        if (add_EM_S) {
3386                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3387                            if (d.actsExpanded()) {
3388                                const double d_0 = d_p[0];
3389                                const double d_1 = d_p[1];
3390                                const double tmp0_0 = d_0 + d_1;
3391                                const double tmp1_1 = d_1*w1;
3392                                const double tmp4_1 = d_0*w1;
3393                                const double tmp0_1 = d_0*w0;
3394                                const double tmp3_1 = d_1*w0;
3395                                const double tmp2_1 = tmp0_0*w2;
3396                                EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
3397                                EM_S[INDEX2(0,2,4)]+=tmp2_1;
3398                                EM_S[INDEX2(2,0,4)]+=tmp2_1;
3399                                EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;
3400                            } else { /* constant data */
3401                                const double d_0 = d_p[0];
3402                                const double tmp1_1 = d_0*w4;
3403                                const double tmp0_1 = d_0*w3;
3404                                EM_S[INDEX2(0,0,4)]+=tmp0_1;
3405                                EM_S[INDEX2(0,2,4)]+=tmp1_1;
3406                                EM_S[INDEX2(2,0,4)]+=tmp1_1;
3407                                EM_S[INDEX2(2,2,4)]+=tmp0_1;
3408                            }
3409                        }
3410                        ///////////////
3411                        // process y //
3412                        ///////////////
3413                        if (add_EM_F) {
3414                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3415                            if (y.actsExpanded()) {
3416                                const double y_0 = y_p[0];
3417                                const double y_1 = y_p[1];
3418                                const double tmp3_1 = w5*y_1;
3419                                const double tmp2_1 = w6*y_0;
3420                                const double tmp0_1 = w6*y_1;
3421                                const double tmp1_1 = w5*y_0;
3422                                EM_F[0]+=tmp0_1 + tmp1_1;
3423                                EM_F[2]+=tmp2_1 + tmp3_1;
3424                            } else { /* constant data */
3425                                const double y_0 = y_p[0];
3426                                const double tmp0_1 = w7*y_0;
3427                                EM_F[0]+=tmp0_1;
3428                                EM_F[2]+=tmp0_1;
3429                            }
3430                        }
3431                        /* GENERATOR SNIP_PDEBC_SINGLE_0 BOTTOM */
3432                        const index_t firstNode=m_N0*k1;
3433                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3434                    }
3435                } // end colouring
3436            }
3437    
3438            if (m_faceOffset[1] > -1) {
3439                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3440    #pragma omp for nowait
3441                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3442                        vector<double> EM_S(4*4, 0);
3443                        vector<double> EM_F(4, 0);
3444                        const index_t e = m_faceOffset[1]+k1;
3445                        /* GENERATOR SNIP_PDEBC_SINGLE_1 TOP */
3446                        ///////////////
3447                        // process d //
3448                        ///////////////
3449                        if (add_EM_S) {
3450                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3451                            if (d.actsExpanded()) {
3452                                const double d_0 = d_p[0];
3453                                const double d_1 = d_p[1];
3454                                const double tmp0_0 = d_0 + d_1;
3455                                const double tmp4_1 = d_1*w1;
3456                                const double tmp1_1 = d_0*w1;
3457                                const double tmp3_1 = d_0*w0;
3458                                const double tmp0_1 = d_1*w0;
3459                                const double tmp2_1 = tmp0_0*w2;
3460                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp1_1;
3461                                EM_S[INDEX2(3,1,4)]+=tmp2_1;
3462                                EM_S[INDEX2(1,3,4)]+=tmp2_1;
3463                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;
3464                            } else { /* constant data */
3465                                const double d_0 = d_p[0];
3466                                const double tmp1_1 = d_0*w4;
3467                                const double tmp0_1 = d_0*w3;
3468                                EM_S[INDEX2(3,3,4)]+=tmp0_1;
3469                                EM_S[INDEX2(3,1,4)]+=tmp1_1;
3470                                EM_S[INDEX2(1,3,4)]+=tmp1_1;
3471                                EM_S[INDEX2(1,1,4)]+=tmp0_1;
3472                            }
3473                        }
3474                        ///////////////
3475                        // process y //
3476                        ///////////////
3477                        if (add_EM_F) {
3478                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3479                            if (y.actsExpanded()) {
3480                                const double y_0 = y_p[0];
3481                                const double y_1 = y_p[1];
3482                                const double tmp3_1 = w5*y_1;
3483                                const double tmp2_1 = w6*y_0;
3484                                const double tmp0_1 = w6*y_1;
3485                                const double tmp1_1 = w5*y_0;
3486                                EM_F[1]+=tmp0_1 + tmp1_1;
3487                                EM_F[3]+=tmp2_1 + tmp3_1;
3488                            } else { /* constant data */
3489                                const double y_0 = y_p[0];
3490                                const double tmp0_1 = w7*y_0;
3491                                EM_F[1]+=tmp0_1;
3492                                EM_F[3]+=tmp0_1;
3493                            }
3494                        }
3495                        /* GENERATOR SNIP_PDEBC_SINGLE_1 BOTTOM */
3496                        const index_t firstNode=m_N0*(k1+1)-2;
3497                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3498                    }
3499                } // end colouring
3500            }
3501    
3502            if (m_faceOffset[2] > -1) {
3503                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3504    #pragma omp for nowait
3505                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3506                        vector<double> EM_S(4*4, 0);
3507                        vector<double> EM_F(4, 0);
3508                        const index_t e = m_faceOffset[2]+k0;
3509                        /* GENERATOR SNIP_PDEBC_SINGLE_2 TOP */
3510                        ///////////////
3511                        // process d //
3512                        ///////////////
3513                        if (add_EM_S) {
3514                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3515                            if (d.actsExpanded()) {
3516                                const double d_0 = d_p[0];
3517                                const double d_1 = d_p[1];
3518                                const double tmp0_0 = d_0 + d_1;
3519                                const double tmp4_1 = d_1*w9;
3520                                const double tmp2_1 = d_0*w9;
3521                                const double tmp0_1 = tmp0_0*w8;
3522                                const double tmp1_1 = d_1*w10;
3523                                const double tmp3_1 = d_0*w10;
3524                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
3525                                EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;
3526                                EM_S[INDEX2(1,0,4)]+=tmp0_1;
3527                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;
3528                            } else { /* constant data */
3529                                const double d_0 = d_p[0];
3530                                const double tmp0_1 = d_0*w11;
3531                                const double tmp1_1 = d_0*w12;
3532                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
3533                                EM_S[INDEX2(0,0,4)]+=tmp1_1;
3534                                EM_S[INDEX2(1,0,4)]+=tmp0_1;
3535                                EM_S[INDEX2(1,1,4)]+=tmp1_1;
3536                            }
3537                        }
3538                        ///////////////
3539                        // process y //
3540                        ///////////////
3541                        if (add_EM_F) {
3542                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3543                            if (y.actsExpanded()) {
3544                                const double y_0 = y_p[0];
3545                                const double y_1 = y_p[1];
3546                                const double tmp2_1 = w13*y_1;
3547                                const double tmp1_1 = w14*y_1;
3548                                const double tmp3_1 = w14*y_0;
3549                                const double tmp0_1 = w13*y_0;
3550                                EM_F[0]+=tmp0_1 + tmp1_1;
3551                                EM_F[1]+=tmp2_1 + tmp3_1;
3552                            } else { /* constant data */
3553                                const double y_0 = y_p[0];
3554                                const double tmp0_1 = w15*y_0;
3555                                EM_F[0]+=tmp0_1;
3556                                EM_F[1]+=tmp0_1;
3557                            }
3558                        }
3559                        /* GENERATOR SNIP_PDEBC_SINGLE_2 BOTTOM */
3560                        const index_t firstNode=k0;
3561                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3562                    }
3563                } // end colouring
3564            }
3565    
3566            if (m_faceOffset[3] > -1) {
3567                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3568    #pragma omp for nowait
3569                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3570                        const index_t e = m_faceOffset[3]+k0;
3571                        vector<double> EM_S(4*4, 0);
3572                        vector<double> EM_F(4, 0);
3573                        /* GENERATOR SNIP_PDEBC_SINGLE_3 TOP */
3574                        ///////////////
3575                        // process d //
3576                        ///////////////
3577                        if (add_EM_S) {
3578                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3579                            if (d.actsExpanded()) {
3580                                const double d_0 = d_p[0];
3581                                const double d_1 = d_p[1];
3582                                const double tmp0_0 = d_0 + d_1;
3583                                const double tmp2_1 = d_1*w9;
3584                                const double tmp4_1 = d_0*w9;
3585                                const double tmp0_1 = tmp0_0*w8;
3586                                const double tmp3_1 = d_1*w10;
3587                                const double tmp1_1 = d_0*w10;
3588                                EM_S[INDEX2(3,2,4)]+=tmp0_1;
3589                                EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
3590                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
3591                                EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;
3592                            } else { /* constant data */
3593                                const double d_0 = d_p[0];
3594                                const double tmp0_1 = d_0*w11;
3595                                const double tmp1_1 = d_0*w12;
3596                                EM_S[INDEX2(3,2,4)]+=tmp0_1;
3597                                EM_S[INDEX2(3,3,4)]+=tmp1_1;
3598                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
3599                                EM_S[INDEX2(2,2,4)]+=tmp1_1;
3600                            }
3601                        }
3602                        ///////////////
3603                        // process y //
3604                        ///////////////
3605                        if (add_EM_F) {
3606                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3607                            if (y.actsExpanded()) {
3608                                const double y_0 = y_p[0];
3609                                const double y_1 = y_p[1];
3610                                const double tmp2_1 = w13*y_1;
3611                                const double tmp1_1 = w14*y_1;
3612                                const double tmp3_1 = w14*y_0;
3613                                const double tmp0_1 = w13*y_0;
3614                                EM_F[2]+=tmp0_1 + tmp1_1;
3615                                EM_F[3]+=tmp2_1 + tmp3_1;
3616                            } else { /* constant data */
3617                                const double y_0 = y_p[0];
3618                                const double tmp0_1 = w15*y_0;
3619                                EM_F[2]+=tmp0_1;
3620                                EM_F[3]+=tmp0_1;
3621                            }
3622                        }
3623                        /* GENERATOR SNIP_PDEBC_SINGLE_3 BOTTOM */
3624                        const index_t firstNode=m_N0*(m_N1-2)+k0;
3625                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3626                    }
3627                } // end colouring
3628            }
3629        } // end of parallel section
3630    }
3631    
3632    //protected
3633    void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
3634          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3635    {
3636        const double h0 = m_l0/m_gNE0;
3637        const double h1 = m_l1/m_gNE1;
3638        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE TOP */
3639        const double w0 = 0.25*h1;
3640        const double w1 = 0.5*h1;
3641        const double w2 = 0.25*h0;
3642        const double w3 = 0.5*h0;
3643        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE BOTTOM */
3644        const bool add_EM_S=!d.isEmpty();
3645        const bool add_EM_F=!y.isEmpty();
3646    #pragma omp parallel
3647        {
3648            if (m_faceOffset[0] > -1) {
3649                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3650    #pragma omp for nowait
3651                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3652                        vector<double> EM_S(4*4, 0);
3653                        vector<double> EM_F(4, 0);
3654                        const index_t e = k1;
3655                        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 TOP */
3656                        ///////////////
3657                        // process d //
3658                        ///////////////
3659                        if (add_EM_S) {
3660                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3661                            const double d_0 = d_p[0];
3662                            const double tmp0_1 = d_0*w0;
3663                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
3664                            EM_S[INDEX2(0,2,4)]+=tmp0_1;
3665                            EM_S[INDEX2(2,0,4)]+=tmp0_1;
3666                            EM_S[INDEX2(2,2,4)]+=tmp0_1;
3667                        }
3668                        ///////////////
3669                        // process y //
3670                        ///////////////
3671                        if (add_EM_F) {
3672                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3673                            const double y_0 = y_p[0];
3674                            const double tmp0_1 = w1*y_0;
3675                            EM_F[0]+=tmp0_1;
3676                            EM_F[2]+=tmp0_1;
3677                        }
3678                        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 BOTTOM */
3679                        const index_t firstNode=m_N0*k1;
3680                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3681                    }
3682                } // end colouring
3683            }
3684    
3685            if (m_faceOffset[1] > -1) {
3686                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3687    #pragma omp for nowait
3688                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3689                        vector<double> EM_S(4*4, 0);
3690                        vector<double> EM_F(4, 0);
3691                        const index_t e = m_faceOffset[1]+k1;
3692                        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 TOP */
3693                        ///////////////
3694                        // process d //
3695                        ///////////////
3696                        if (add_EM_S) {
3697                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3698                            const double d_0 = d_p[0];
3699                            const double tmp0_1 = d_0*w0;
3700                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
3701                            EM_S[INDEX2(3,1,4)]+=tmp0_1;
3702                            EM_S[INDEX2(1,3,4)]+=tmp0_1;
3703                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
3704                        }
3705                        ///////////////
3706                        // process y //
3707                        ///////////////
3708                        if (add_EM_F) {
3709                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3710                            const double y_0 = y_p[0];
3711                            const double tmp0_1 = w1*y_0;
3712                            EM_F[1]+=tmp0_1;
3713                            EM_F[3]+=tmp0_1;
3714                        }
3715                        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 BOTTOM */
3716                        const index_t firstNode=m_N0*(k1+1)-2;
3717                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3718                    }
3719                } // end colouring
3720            }
3721    
3722            if (m_faceOffset[2] > -1) {
3723                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3724    #pragma omp for nowait
3725                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3726                        vector<double> EM_S(4*4, 0);
3727                        vector<double> EM_F(4, 0);
3728                        const index_t e = m_faceOffset[2]+k0;
3729                        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 TOP */
3730                        ///////////////
3731                        // process d //
3732                        ///////////////
3733                        if (add_EM_S) {
3734                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3735                            const double d_0 = d_p[0];
3736                            const double tmp0_1 = d_0*w2;
3737                            EM_S[INDEX2(0,1,4)]+=tmp0_1;
3738                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
3739                            EM_S[INDEX2(1,0,4)]+=tmp0_1;
3740                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
3741                        }
3742                        ///////////////
3743                        // process y //
3744                        ///////////////
3745                        if (add_EM_F) {
3746                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3747                            const double y_0 = y_p[0];
3748                            const double tmp0_1 = w3*y_0;
3749                            EM_F[0]+=tmp0_1;
3750                            EM_F[1]+=tmp0_1;
3751                        }
3752                        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 BOTTOM */
3753                        const index_t firstNode=k0;
3754                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3755                    }
3756                } // end colouring
3757            }
3758    
3759            if (m_faceOffset[3] > -1) {
3760                for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3761    #pragma omp for nowait
3762                    for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3763                        vector<double> EM_S(4*4, 0);
3764                        vector<double> EM_F(4, 0);
3765                        const index_t e = m_faceOffset[3]+k0;
3766                        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 TOP */
3767                        ///////////////
3768                        // process d //
3769                        ///////////////
3770                        if (add_EM_S) {
3771                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3772                            const double d_0 = d_p[0];
3773                            const double tmp0_1 = d_0*w2;
3774                            EM_S[INDEX2(3,2,4)]+=tmp0_1;
3775                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
3776                            EM_S[INDEX2(2,3,4)]+=tmp0_1;
3777                            EM_S[INDEX2(2,2,4)]+=tmp0_1;
3778                        }
3779                        ///////////////
3780                        // process y //
3781                        ///////////////
3782                        if (add_EM_F) {
3783                            const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3784                            const double y_0 = y_p[0];
3785                            const double tmp0_1 = w3*y_0;
3786                            EM_F[2]+=tmp0_1;
3787                            EM_F[3]+=tmp0_1;
3788                        }
3789                        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 BOTTOM */
3790                        const index_t firstNode=m_N0*(m_N1-2)+k0;
3791                        addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
3792                    }
3793                } // end colouring
3794            }
3795        } // end of parallel section
3796    }
3797    
3798    //protected
3799    void Rectangle::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
3800          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3801    {
3802        const double h0 = m_l0/m_gNE0;
3803        const double h1 = m_l1/m_gNE1;
3804        dim_t numEq, numComp;
3805        if (!mat) {
3806            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
3807        } else {
3808            numEq=mat->logical_row_block_size;
3809            numComp=mat->logical_col_block_size;
3810        }
3811        /* GENERATOR SNIP_PDEBC_SYSTEM_PRE TOP */
3812        const double w0 = 0.31100423396407310779*h1;
3813        const double w1 = 0.022329099369260225539*h1;
3814        const double w10 = 0.022329099369260225539*h0;
3815        const double w11 = 0.16666666666666666667*h0;
3816        const double w12 = 0.33333333333333333333*h0;
3817        const double w13 = 0.39433756729740644113*h0;
3818        const double w14 = 0.10566243270259355887*h0;
3819        const double w15 = 0.5*h0;
3820        const double w2 = 0.083333333333333333333*h1;
3821        const double w3 = 0.33333333333333333333*h1;
3822        const double w4 = 0.16666666666666666667*h1;
3823        const double w5 = 0.39433756729740644113*h1;
3824        const double w6 = 0.10566243270259355887*h1;
3825        const double w7 = 0.5*h1;
3826        const double w8 = 0.083333333333333333333*h0;
3827        const double w9 = 0.31100423396407310779*h0;
3828        /* GENERATOR SNIP_PDEBC_SYSTEM_PRE BOTTOM */
3829        const bool add_EM_S=!d.isEmpty();
3830        const bool add_EM_F=!y.isEmpty();
3831    #pragma omp parallel
3832        {
3833            if (m_faceOffset[0] > -1) {
3834                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3835    #pragma omp for nowait
3836                    for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3837                        vector<double> EM_S(4*4*numEq*numComp, 0);
3838                        vector<double> EM_F(4*numEq, 0);
3839                        const index_t e = k1;
3840                        /* GENERATOR SNIP_PDEBC_SYSTEM_0 TOP */
3841                        ///////////////
3842                        // process d //
3843                        ///////////////
3844                        if (add_EM_S) {
3845                            const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3846                            if (d.actsExpanded()) {
3847                                for (index_t k=0; k<numEq; k++) {
3848                                    for (index_t m=0; m<numComp; m++) {
3849                                        const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];
3850                                        const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];
3851                                        const double tmp0_0 = d_0 + d_1;
3852                                        const double tmp1_1 = d_1*w1;
3853                                        const double tmp4_1 = d_0*w1;
3854                                        const double tmp0_1 = d_0*w0;
3855                                        const double tmp3_1 = d_1*w0;
3856                                        const double tmp2_1 = tmp0_0*w2;