/[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 3768 by caltinay, Fri Jan 13 07:17:16 2012 UTC revision 3769 by caltinay, Tue Jan 17 04:09:55 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 643  void Rectangle::assembleGradient(escript Line 630  void Rectangle::assembleGradient(escript
630  #pragma omp parallel for  #pragma omp parallel for
631          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
632              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
633                  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));
634                  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));
635                  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));
636                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
637                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
638                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
639                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      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 652  void Rectangle::assembleGradient(escript
652  #pragma omp parallel for  #pragma omp parallel for
653          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
654              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
655                  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));
656                  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));
657                  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));
658                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
659                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
660                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
661                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                      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 670  void Rectangle::assembleGradient(escript
670              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
671  #pragma omp for nowait  #pragma omp for nowait
672                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
673                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
674                      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));
675                      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));
676                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
677                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
678                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
679                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          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 686  void Rectangle::assembleGradient(escript
686              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
687  #pragma omp for nowait  #pragma omp for nowait
688                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
689                      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));
690                      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));
691                      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));
692                      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));
693                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
694                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
695                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          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 702  void Rectangle::assembleGradient(escript
702              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
703  #pragma omp for nowait  #pragma omp for nowait
704                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
705                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
706                      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));
707                      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));
708                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
709                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
710                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
711                          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 718  void Rectangle::assembleGradient(escript
718              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
719  #pragma omp for nowait  #pragma omp for nowait
720                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
721                      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));
722                      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));
723                      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));
724                      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));
725                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
726                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
727                          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 739  void Rectangle::assembleGradient(escript
739              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
740  #pragma omp for nowait  #pragma omp for nowait
741                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
742                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
743                      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));
744                      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));
745                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
746                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
747                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
748                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          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 753  void Rectangle::assembleGradient(escript
753              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
754  #pragma omp for nowait  #pragma omp for nowait
755                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
756                      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));
757                      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));
758                      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));
759                      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));
760                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
761                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
762                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          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 767  void Rectangle::assembleGradient(escript
767              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
768  #pragma omp for nowait  #pragma omp for nowait
769                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
770                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
771                      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));
772                      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));
773                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
774                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
775                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
776                          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 781  void Rectangle::assembleGradient(escript
781              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
782  #pragma omp for nowait  #pragma omp for nowait
783                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
784                      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));
785                      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));
786                      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));
787                      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));
788                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
789                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
790                          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 814  void Rectangle::assembleIntegrate(vector
814                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
815                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
816                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
817                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
818                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
819                          const register double f2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
820                          const register double f3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
821                          int_local[i]+=(f0+f1+f2+f3)*w;                          int_local[i]+=(f0+f1+f2+f3)*w;
822                      }  /* end of component loop i */                      }  /* end of component loop i */
823                  } /* end of k0 loop */                  } /* end of k0 loop */
# Line 870  void Rectangle::assembleIntegrate(vector Line 857  void Rectangle::assembleIntegrate(vector
857                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
858                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
859                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
860                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
861                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
862                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
863                      }  /* end of component loop i */                      }  /* end of component loop i */
864                  } /* end of k1 loop */                  } /* end of k1 loop */
# Line 882  void Rectangle::assembleIntegrate(vector Line 869  void Rectangle::assembleIntegrate(vector
869                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
870                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
871                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
872                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
873                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
874                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
875                      }  /* end of component loop i */                      }  /* end of component loop i */
876                  } /* end of k1 loop */                  } /* end of k1 loop */
# Line 894  void Rectangle::assembleIntegrate(vector Line 881  void Rectangle::assembleIntegrate(vector
881                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
882                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
883                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
884                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
885                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
886                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
887                      }  /* end of component loop i */                      }  /* end of component loop i */
888                  } /* end of k0 loop */                  } /* end of k0 loop */
# Line 906  void Rectangle::assembleIntegrate(vector Line 893  void Rectangle::assembleIntegrate(vector
893                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
894                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
895                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
896                          const register double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
897                          const register double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
898                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
899                      }  /* end of component loop i */                      }  /* end of component loop i */
900                  } /* end of k0 loop */                  } /* end of k0 loop */
# Line 1296  void Rectangle::interpolateNodesOnElemen Line 1283  void Rectangle::interpolateNodesOnElemen
1283  #pragma omp parallel for  #pragma omp parallel for
1284          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1285              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
1286                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1287                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1288                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1289                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1290                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1291                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1292                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
# Line 1314  void Rectangle::interpolateNodesOnElemen Line 1301  void Rectangle::interpolateNodesOnElemen
1301  #pragma omp parallel for  #pragma omp parallel for
1302          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1303              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
1304                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1305                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1306                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1307                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1308                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1309                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
1310                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);
# Line 1343  void Rectangle::interpolateNodesOnFaces( Line 1330  void Rectangle::interpolateNodesOnFaces(
1330              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1331  #pragma omp for nowait  #pragma omp for nowait
1332                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1333                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1334                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1335                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1336                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1337                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
# Line 1354  void Rectangle::interpolateNodesOnFaces( Line 1341  void Rectangle::interpolateNodesOnFaces(
1341              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1342  #pragma omp for nowait  #pragma omp for nowait
1343                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1344                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1345                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1346                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1347                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1348                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
# Line 1365  void Rectangle::interpolateNodesOnFaces( Line 1352  void Rectangle::interpolateNodesOnFaces(
1352              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1353  #pragma omp for nowait  #pragma omp for nowait
1354                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1355                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1356                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1357                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1358                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1359                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
# Line 1376  void Rectangle::interpolateNodesOnFaces( Line 1363  void Rectangle::interpolateNodesOnFaces(
1363              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1364  #pragma omp for nowait  #pragma omp for nowait
1365                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1366                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1367                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1368                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1369                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1370                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
# Line 1394  void Rectangle::interpolateNodesOnFaces( Line 1381  void Rectangle::interpolateNodesOnFaces(
1381              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1382  #pragma omp for nowait  #pragma omp for nowait
1383                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1384                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1385                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1386                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1387                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1388                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;
# Line 1406  void Rectangle::interpolateNodesOnFaces( Line 1393  void Rectangle::interpolateNodesOnFaces(
1393              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1394  #pragma omp for nowait  #pragma omp for nowait
1395                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1396                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1397                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1398                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1399                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1400                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;
# Line 1418  void Rectangle::interpolateNodesOnFaces( Line 1405  void Rectangle::interpolateNodesOnFaces(
1405              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1406  #pragma omp for nowait  #pragma omp for nowait
1407                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1408                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1409                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1410                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1411                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1412                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;
# Line 1430  void Rectangle::interpolateNodesOnFaces( Line 1417  void Rectangle::interpolateNodesOnFaces(
1417              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1418  #pragma omp for nowait  #pragma omp for nowait
1419                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1420                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1421                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1422                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1423                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1424                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;
# Line 1447  void Rectangle::interpolateNodesOnFaces( Line 1434  void Rectangle::interpolateNodesOnFaces(
1434  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1435          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1436          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
1437          const escript::Data& X, const escript::Data& Y,          const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
1438  {  {
1439      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1440      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 1471  void Rectangle::assemblePDESingle(Paso_S Line 1457  void Rectangle::assemblePDESingle(Paso_S
1457      const double w16 = -0.01116454968463011277*h0/h1;      const double w16 = -0.01116454968463011277*h0/h1;
1458      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*h0/h1;
1459      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*h1/h0;
1460      const double w19 = 0.25000000000000000000;      const double w19 = 0.25;
1461      const double w20 = -0.25000000000000000000;      const double w20 = -0.25;
1462      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*h0/h1;
1463      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*h1/h0;
1464      const double w23 = -0.16666666666666666667*h0/h1;      const double w23 = -0.16666666666666666667*h0/h1;
# Line 1548  void Rectangle::assemblePDESingle(Paso_S Line 1534  void Rectangle::assemblePDESingle(Paso_S
1534                          add_EM_S=true;                          add_EM_S=true;
1535                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1536                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
1537                              const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];                              const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1538                              const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];                              const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1539                              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)];
1540                              const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];                              const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1541                              const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];                              const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1542                              const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];                              const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1543                              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)];
1544                              const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];                              const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1545                              const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];                              const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1546                              const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];                              const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1547                              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)];
1548                              const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];                              const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1549                              const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];                              const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1550                              const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];                              const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1551                              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)];
1552                              const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];                              const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1553                              const register double tmp0_0 = A_01_0 + A_01_3;                              const double tmp0_0 = A_01_0 + A_01_3;
1554                              const register double tmp1_0 = A_00_0 + A_00_1;                              const double tmp1_0 = A_00_0 + A_00_1;
1555                              const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                              const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1556                              const register double tmp3_0 = A_00_2 + A_00_3;                              const double tmp3_0 = A_00_2 + A_00_3;
1557                              const register double tmp4_0 = A_10_1 + A_10_2;                              const double tmp4_0 = A_10_1 + A_10_2;
1558                              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;
1559                              const register double tmp6_0 = A_01_3 + A_10_0;                              const double tmp6_0 = A_01_3 + A_10_0;
1560                              const register double tmp7_0 = A_01_0 + A_10_3;                              const double tmp7_0 = A_01_0 + A_10_3;
1561                              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;
1562                              const register double tmp9_0 = A_01_0 + A_10_0;                              const double tmp9_0 = A_01_0 + A_10_0;
1563                              const register double tmp12_0 = A_11_0 + A_11_2;                              const double tmp12_0 = A_11_0 + A_11_2;
1564                              const register double tmp10_0 = A_01_3 + A_10_3;                              const double tmp10_0 = A_01_3 + A_10_3;
1565                              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;
1566                              const register double tmp13_0 = A_01_2 + A_10_1;                              const double tmp13_0 = A_01_2 + A_10_1;
1567                              const register double tmp11_0 = A_11_1 + A_11_3;                              const double tmp11_0 = A_11_1 + A_11_3;
1568                              const register double tmp18_0 = A_01_1 + A_10_1;                              const double tmp18_0 = A_01_1 + A_10_1;
1569                              const register double tmp15_0 = A_01_1 + A_10_2;                              const double tmp15_0 = A_01_1 + A_10_2;
1570                              const register double tmp16_0 = A_10_0 + A_10_3;                              const double tmp16_0 = A_10_0 + A_10_3;
1571                              const register double tmp17_0 = A_01_1 + A_01_2;                              const double tmp17_0 = A_01_1 + A_01_2;
1572                              const register double tmp19_0 = A_01_2 + A_10_2;                              const double tmp19_0 = A_01_2 + A_10_2;
1573                              const register double tmp0_1 = A_10_3*w8;                              const double tmp0_1 = A_10_3*w8;
1574                              const register double tmp1_1 = tmp0_0*w1;                              const double tmp1_1 = tmp0_0*w1;
1575                              const register double tmp2_1 = A_01_1*w4;                              const double tmp2_1 = A_01_1*w4;
1576                              const register double tmp3_1 = tmp1_0*w0;                              const double tmp3_1 = tmp1_0*w0;
1577                              const register double tmp4_1 = A_01_2*w7;                              const double tmp4_1 = A_01_2*w7;
1578                              const register double tmp5_1 = tmp2_0*w3;                              const double tmp5_1 = tmp2_0*w3;
1579                              const register double tmp6_1 = tmp3_0*w6;                              const double tmp6_1 = tmp3_0*w6;
1580                              const register double tmp7_1 = A_10_0*w2;                              const double tmp7_1 = A_10_0*w2;
1581                              const register double tmp8_1 = tmp4_0*w5;                              const double tmp8_1 = tmp4_0*w5;
1582                              const register double tmp9_1 = tmp2_0*w10;                              const double tmp9_1 = tmp2_0*w10;
1583                              const register double tmp14_1 = A_10_0*w8;                              const double tmp14_1 = A_10_0*w8;
1584                              const register double tmp23_1 = tmp3_0*w14;                              const double tmp23_1 = tmp3_0*w14;
1585                              const register double tmp35_1 = A_01_0*w8;                              const double tmp35_1 = A_01_0*w8;
1586                              const register double tmp54_1 = tmp13_0*w8;                              const double tmp54_1 = tmp13_0*w8;
1587                              const register double tmp20_1 = tmp9_0*w4;                              const double tmp20_1 = tmp9_0*w4;
1588                              const register double tmp25_1 = tmp12_0*w12;                              const double tmp25_1 = tmp12_0*w12;
1589                              const register double tmp44_1 = tmp7_0*w7;                              const double tmp44_1 = tmp7_0*w7;
1590                              const register double tmp26_1 = tmp10_0*w4;                              const double tmp26_1 = tmp10_0*w4;
1591                              const register double tmp52_1 = tmp18_0*w8;                              const double tmp52_1 = tmp18_0*w8;
1592                              const register double tmp48_1 = A_10_1*w7;                              const double tmp48_1 = A_10_1*w7;
1593                              const register double tmp46_1 = A_01_3*w8;                              const double tmp46_1 = A_01_3*w8;
1594                              const register double tmp50_1 = A_01_0*w2;                              const double tmp50_1 = A_01_0*w2;
1595                              const register double tmp56_1 = tmp19_0*w8;                              const double tmp56_1 = tmp19_0*w8;
1596                              const register double tmp19_1 = A_10_3*w2;                              const double tmp19_1 = A_10_3*w2;
1597                              const register double tmp47_1 = A_10_2*w4;                              const double tmp47_1 = A_10_2*w4;
1598                              const register double tmp16_1 = tmp3_0*w0;                              const double tmp16_1 = tmp3_0*w0;
1599                              const register double tmp18_1 = tmp1_0*w6;                              const double tmp18_1 = tmp1_0*w6;
1600                              const register double tmp31_1 = tmp11_0*w12;                              const double tmp31_1 = tmp11_0*w12;
1601                              const register double tmp55_1 = tmp15_0*w2;                              const double tmp55_1 = tmp15_0*w2;
1602                              const register double tmp39_1 = A_10_2*w7;                              const double tmp39_1 = A_10_2*w7;
1603                              const register double tmp11_1 = tmp6_0*w7;                              const double tmp11_1 = tmp6_0*w7;
1604                              const register double tmp40_1 = tmp11_0*w17;                              const double tmp40_1 = tmp11_0*w17;
1605                              const register double tmp34_1 = tmp15_0*w8;                              const double tmp34_1 = tmp15_0*w8;
1606                              const register double tmp33_1 = tmp14_0*w5;                              const double tmp33_1 = tmp14_0*w5;
1607                              const register double tmp24_1 = tmp11_0*w13;                              const double tmp24_1 = tmp11_0*w13;
1608                              const register double tmp43_1 = tmp17_0*w5;                              const double tmp43_1 = tmp17_0*w5;
1609                              const register double tmp15_1 = A_01_2*w4;                              const double tmp15_1 = A_01_2*w4;
1610                              const register double tmp53_1 = tmp19_0*w2;                              const double tmp53_1 = tmp19_0*w2;
1611                              const register double tmp27_1 = tmp3_0*w11;                              const double tmp27_1 = tmp3_0*w11;
1612                              const register double tmp32_1 = tmp13_0*w2;                              const double tmp32_1 = tmp13_0*w2;
1613                              const register double tmp10_1 = tmp5_0*w9;                              const double tmp10_1 = tmp5_0*w9;
1614                              const register double tmp37_1 = A_10_1*w4;                              const double tmp37_1 = A_10_1*w4;
1615                              const register double tmp38_1 = tmp5_0*w15;                              const double tmp38_1 = tmp5_0*w15;
1616                              const register double tmp17_1 = A_01_1*w7;                              const double tmp17_1 = A_01_1*w7;
1617                              const register double tmp12_1 = tmp7_0*w4;                              const double tmp12_1 = tmp7_0*w4;
1618                              const register double tmp22_1 = tmp10_0*w7;                              const double tmp22_1 = tmp10_0*w7;
1619                              const register double tmp57_1 = tmp18_0*w2;                              const double tmp57_1 = tmp18_0*w2;
1620                              const register double tmp28_1 = tmp9_0*w7;                              const double tmp28_1 = tmp9_0*w7;
1621                              const register double tmp29_1 = tmp1_0*w14;                              const double tmp29_1 = tmp1_0*w14;
1622                              const register double tmp51_1 = tmp11_0*w16;                              const double tmp51_1 = tmp11_0*w16;
1623                              const register double tmp42_1 = tmp12_0*w16;                              const double tmp42_1 = tmp12_0*w16;
1624                              const register double tmp49_1 = tmp12_0*w17;                              const double tmp49_1 = tmp12_0*w17;
1625                              const register double tmp21_1 = tmp1_0*w11;                              const double tmp21_1 = tmp1_0*w11;
1626                              const register double tmp45_1 = tmp6_0*w4;                              const double tmp45_1 = tmp6_0*w4;
1627                              const register double tmp13_1 = tmp8_0*w1;                              const double tmp13_1 = tmp8_0*w1;
1628                              const register double tmp36_1 = tmp16_0*w1;                              const double tmp36_1 = tmp16_0*w1;
1629                              const register double tmp41_1 = A_01_3*w2;                              const double tmp41_1 = A_01_3*w2;
1630                              const register double tmp30_1 = tmp12_0*w13;                              const double tmp30_1 = tmp12_0*w13;
1631                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1632                              EM_S[INDEX2(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;
1633                              EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;                              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 1645  void Rectangle::assemblePDESingle(Paso_S
1645                              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;
1646                              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;
1647                          } else { // constant data                          } else { // constant data
1648                              const register double A_00 = A_p[INDEX2(0,0,2)];                              const double A_00 = A_p[INDEX2(0,0,2)];
1649                              const register double A_10 = A_p[INDEX2(1,0,2)];                              const double A_10 = A_p[INDEX2(1,0,2)];
1650                              const register double A_01 = A_p[INDEX2(0,1,2)];                              const double A_01 = A_p[INDEX2(0,1,2)];
1651                              const register double A_11 = A_p[INDEX2(1,1,2)];                              const double A_11 = A_p[INDEX2(1,1,2)];
1652                              const register double tmp0_0 = A_01 + A_10;                              const double tmp0_0 = A_01 + A_10;
1653                              const register double tmp0_1 = A_00*w18;                              const double tmp0_1 = A_00*w18;
1654                              const register double tmp1_1 = A_01*w19;                              const double tmp1_1 = A_01*w19;
1655                              const register double tmp2_1 = A_10*w20;                              const double tmp2_1 = A_10*w20;
1656                              const register double tmp3_1 = A_11*w21;                              const double tmp3_1 = A_11*w21;
1657                              const register double tmp4_1 = A_00*w22;                              const double tmp4_1 = A_00*w22;
1658                              const register double tmp5_1 = tmp0_0*w19;                              const double tmp5_1 = tmp0_0*w19;
1659                              const register double tmp6_1 = A_11*w23;                              const double tmp6_1 = A_11*w23;
1660                              const register double tmp7_1 = A_11*w25;                              const double tmp7_1 = A_11*w25;
1661                              const register double tmp8_1 = A_00*w24;                              const double tmp8_1 = A_00*w24;
1662                              const register double tmp9_1 = tmp0_0*w20;                              const double tmp9_1 = tmp0_0*w20;
1663                              const register double tmp10_1 = A_01*w20;                              const double tmp10_1 = A_01*w20;
1664                              const register double tmp11_1 = A_11*w27;                              const double tmp11_1 = A_11*w27;
1665                              const register double tmp12_1 = A_00*w26;                              const double tmp12_1 = A_00*w26;
1666                              const register double tmp13_1 = A_10*w19;                              const double tmp13_1 = A_10*w19;
1667                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1668                              EM_S[INDEX2(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;
1669                              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 1689  void Rectangle::assemblePDESingle(Paso_S
1689                          add_EM_S=true;                          add_EM_S=true;
1690                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1691                          if (B.actsExpanded()) {                          if (B.actsExpanded()) {
1692                              const register double B_0_0 = B_p[INDEX2(0,0,2)];                              const double B_0_0 = B_p[INDEX2(0,0,2)];
1693                              const register double B_1_0 = B_p[INDEX2(1,0,2)];                              const double B_1_0 = B_p[INDEX2(1,0,2)];
1694                              const register double B_0_1 = B_p[INDEX2(0,1,2)];                              const double B_0_1 = B_p[INDEX2(0,1,2)];
1695                              const register double B_1_1 = B_p[INDEX2(1,1,2)];                              const double B_1_1 = B_p[INDEX2(1,1,2)];
1696                              const register double B_0_2 = B_p[INDEX2(0,2,2)];                              const double B_0_2 = B_p[INDEX2(0,2,2)];
1697                              const register double B_1_2 = B_p[INDEX2(1,2,2)];                              const double B_1_2 = B_p[INDEX2(1,2,2)];
1698                              const register double B_0_3 = B_p[INDEX2(0,3,2)];                              const double B_0_3 = B_p[INDEX2(0,3,2)];
1699                              const register double B_1_3 = B_p[INDEX2(1,3,2)];                              const double B_1_3 = B_p[INDEX2(1,3,2)];
1700                              const register double tmp0_0 = B_1_0 + B_1_1;                              const double tmp0_0 = B_1_0 + B_1_1;
1701                              const register double tmp1_0 = B_1_2 + B_1_3;                              const double tmp1_0 = B_1_2 + B_1_3;
1702                              const register double tmp2_0 = B_0_1 + B_0_3;                              const double tmp2_0 = B_0_1 + B_0_3;
1703                              const register double tmp3_0 = B_0_0 + B_0_2;                              const double tmp3_0 = B_0_0 + B_0_2;
1704                              const register double tmp63_1 = B_1_1*w42;                              const double tmp63_1 = B_1_1*w42;
1705                              const register double tmp79_1 = B_1_1*w40;                              const double tmp79_1 = B_1_1*w40;
1706                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
1707                              const register double tmp8_1 = tmp0_0*w32;                              const double tmp8_1 = tmp0_0*w32;
1708                              const register double tmp71_1 = B_0_1*w34;                              const double tmp71_1 = B_0_1*w34;
1709                              const register double tmp19_1 = B_0_3*w31;                              const double tmp19_1 = B_0_3*w31;
1710                              const register double tmp15_1 = B_0_3*w34;                              const double tmp15_1 = B_0_3*w34;
1711                              const register double tmp9_1 = tmp3_0*w34;                              const double tmp9_1 = tmp3_0*w34;
1712                              const register double tmp35_1 = B_1_0*w36;                              const double tmp35_1 = B_1_0*w36;
1713                              const register double tmp66_1 = B_0_3*w28;                              const double tmp66_1 = B_0_3*w28;
1714                              const register double tmp28_1 = B_1_0*w42;                              const double tmp28_1 = B_1_0*w42;
1715                              const register double tmp22_1 = B_1_0*w40;                              const double tmp22_1 = B_1_0*w40;
1716                              const register double tmp16_1 = B_1_2*w29;                              const double tmp16_1 = B_1_2*w29;
1717                              const register double tmp6_1 = tmp2_0*w35;                              const double tmp6_1 = tmp2_0*w35;
1718                              const register double tmp55_1 = B_1_3*w40;                              const double tmp55_1 = B_1_3*w40;
1719                              const register double tmp50_1 = B_1_3*w42;                              const double tmp50_1 = B_1_3*w42;
1720                              const register double tmp7_1 = tmp1_0*w29;                              const double tmp7_1 = tmp1_0*w29;
1721                              const register double tmp1_1 = tmp1_0*w32;                              const double tmp1_1 = tmp1_0*w32;
1722                              const register double tmp57_1 = B_0_3*w30;                              const double tmp57_1 = B_0_3*w30;
1723                              const register double tmp18_1 = B_1_1*w32;                              const double tmp18_1 = B_1_1*w32;
1724                              const register double tmp53_1 = B_1_0*w41;                              const double tmp53_1 = B_1_0*w41;
1725                              const register double tmp61_1 = B_1_3*w36;                              const double tmp61_1 = B_1_3*w36;
1726                              const register double tmp27_1 = B_0_3*w38;                              const double tmp27_1 = B_0_3*w38;
1727                              const register double tmp64_1 = B_0_2*w30;                              const double tmp64_1 = B_0_2*w30;
1728                              const register double tmp76_1 = B_0_1*w38;                              const double tmp76_1 = B_0_1*w38;
1729                              const register double tmp39_1 = tmp2_0*w34;                              const double tmp39_1 = tmp2_0*w34;
1730                              const register double tmp62_1 = B_0_1*w31;                              const double tmp62_1 = B_0_1*w31;
1731                              const register double tmp56_1 = B_0_0*w31;                              const double tmp56_1 = B_0_0*w31;
1732                              const register double tmp49_1 = B_1_1*w36;                              const double tmp49_1 = B_1_1*w36;
1733                              const register double tmp2_1 = B_0_2*w31;                              const double tmp2_1 = B_0_2*w31;
1734                              const register double tmp23_1 = B_0_2*w33;                              const double tmp23_1 = B_0_2*w33;
1735                              const register double tmp38_1 = B_1_1*w43;                              const double tmp38_1 = B_1_1*w43;
1736                              const register double tmp74_1 = B_1_2*w41;                              const double tmp74_1 = B_1_2*w41;
1737                              const register double tmp43_1 = B_1_1*w41;                              const double tmp43_1 = B_1_1*w41;
1738                              const register double tmp58_1 = B_0_2*w28;                              const double tmp58_1 = B_0_2*w28;
1739                              const register double tmp67_1 = B_0_0*w33;                              const double tmp67_1 = B_0_0*w33;
1740                              const register double tmp33_1 = tmp0_0*w39;                              const double tmp33_1 = tmp0_0*w39;
1741                              const register double tmp4_1 = B_0_0*w28;                              const double tmp4_1 = B_0_0*w28;
1742                              const register double tmp20_1 = B_0_0*w30;                              const double tmp20_1 = B_0_0*w30;
1743                              const register double tmp13_1 = B_0_2*w38;                              const double tmp13_1 = B_0_2*w38;
1744                              const register double tmp65_1 = B_1_2*w43;                              const double tmp65_1 = B_1_2*w43;
1745                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
1746                              const register double tmp41_1 = tmp3_0*w33;                              const double tmp41_1 = tmp3_0*w33;
1747                              const register double tmp73_1 = B_0_2*w37;                              const double tmp73_1 = B_0_2*w37;
1748                              const register double tmp69_1 = B_0_0*w38;                              const double tmp69_1 = B_0_0*w38;
1749                              const register double tmp48_1 = B_1_2*w39;                              const double tmp48_1 = B_1_2*w39;
1750                              const register double tmp59_1 = B_0_1*w33;                              const double tmp59_1 = B_0_1*w33;
1751                              const register double tmp17_1 = B_1_3*w41;                              const double tmp17_1 = B_1_3*w41;
1752                              const register double tmp5_1 = B_0_3*w33;                              const double tmp5_1 = B_0_3*w33;
1753                              const register double tmp3_1 = B_0_1*w30;                              const double tmp3_1 = B_0_1*w30;
1754                              const register double tmp21_1 = B_0_1*w28;                              const double tmp21_1 = B_0_1*w28;
1755                              const register double tmp42_1 = B_1_0*w29;                              const double tmp42_1 = B_1_0*w29;
1756                              const register double tmp54_1 = B_1_2*w32;                              const double tmp54_1 = B_1_2*w32;
1757                              const register double tmp60_1 = B_1_0*w39;                              const double tmp60_1 = B_1_0*w39;
1758                              const register double tmp32_1 = tmp1_0*w36;                              const double tmp32_1 = tmp1_0*w36;
1759                              const register double tmp10_1 = B_0_1*w37;                              const double tmp10_1 = B_0_1*w37;
1760                              const register double tmp14_1 = B_0_0*w35;                              const double tmp14_1 = B_0_0*w35;
1761                              const register double tmp29_1 = B_0_1*w35;                              const double tmp29_1 = B_0_1*w35;
1762                              const register double tmp26_1 = B_1_2*w36;                              const double tmp26_1 = B_1_2*w36;
1763                              const register double tmp30_1 = B_1_3*w43;                              const double tmp30_1 = B_1_3*w43;
1764                              const register double tmp70_1 = B_0_2*w35;                              const double tmp70_1 = B_0_2*w35;
1765                              const register double tmp34_1 = B_1_3*w39;                              const double tmp34_1 = B_1_3*w39;
1766                              const register double tmp51_1 = B_1_0*w43;                              const double tmp51_1 = B_1_0*w43;
1767                              const register double tmp31_1 = B_0_2*w34;                              const double tmp31_1 = B_0_2*w34;
1768                              const register double tmp45_1 = tmp3_0*w28;                              const double tmp45_1 = tmp3_0*w28;
1769                              const register double tmp11_1 = tmp1_0*w39;                              const double tmp11_1 = tmp1_0*w39;
1770                              const register double tmp52_1 = B_1_1*w29;                              const double tmp52_1 = B_1_1*w29;
1771                              const register double tmp44_1 = B_1_3*w32;                              const double tmp44_1 = B_1_3*w32;
1772                              const register double tmp25_1 = B_1_1*w39;                              const double tmp25_1 = B_1_1*w39;
1773                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
1774                              const register double tmp72_1 = B_1_3*w29;                              const double tmp72_1 = B_1_3*w29;
1775                              const register double tmp40_1 = tmp2_0*w28;                              const double tmp40_1 = tmp2_0*w28;
1776                              const register double tmp46_1 = B_1_2*w40;                              const double tmp46_1 = B_1_2*w40;
1777                              const register double tmp36_1 = B_1_2*w42;                              const double tmp36_1 = B_1_2*w42;
1778                              const register double tmp24_1 = B_0_0*w37;                              const double tmp24_1 = B_0_0*w37;
1779                              const register double tmp77_1 = B_0_3*w35;                              const double tmp77_1 = B_0_3*w35;
1780                              const register double tmp68_1 = B_0_3*w37;                              const double tmp68_1 = B_0_3*w37;
1781                              const register double tmp78_1 = B_0_0*w34;                              const double tmp78_1 = B_0_0*w34;
1782                              const register double tmp12_1 = tmp0_0*w36;                              const double tmp12_1 = tmp0_0*w36;
1783                              const register double tmp75_1 = B_1_0*w32;                              const double tmp75_1 = B_1_0*w32;
1784                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1785                              EM_S[INDEX2(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;
1786                              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 1798  void Rectangle::assemblePDESingle(Paso_S
1798                              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;
1799                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1800                          } else { // constant data                          } else { // constant data
1801                              const register double B_0 = B_p[0];                              const double B_0 = B_p[0];
1802                              const register double B_1 = B_p[1];                              const double B_1 = B_p[1];
1803                              const register double tmp0_1 = B_0*w44;                              const double tmp0_1 = B_0*w44;
1804                              const register double tmp1_1 = B_1*w45;                              const double tmp1_1 = B_1*w45;
1805                              const register double tmp2_1 = B_0*w46;                              const double tmp2_1 = B_0*w46;
1806                              const register double tmp3_1 = B_0*w47;                              const double tmp3_1 = B_0*w47;
1807                              const register double tmp4_1 = B_1*w48;                              const double tmp4_1 = B_1*w48;
1808                              const register double tmp5_1 = B_1*w49;                              const double tmp5_1 = B_1*w49;
1809                              const register double tmp6_1 = B_1*w50;                              const double tmp6_1 = B_1*w50;
1810                              const register double tmp7_1 = B_0*w51;                              const double tmp7_1 = B_0*w51;
1811                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1812                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1813                              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 1833  void Rectangle::assemblePDESingle(Paso_S
1833                          add_EM_S=true;                          add_EM_S=true;
1834                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
1835                          if (C.actsExpanded()) {                          if (C.actsExpanded()) {
1836                              const register double C_0_0 = C_p[INDEX2(0,0,2)];                              const double C_0_0 = C_p[INDEX2(0,0,2)];
1837                              const register double C_1_0 = C_p[INDEX2(1,0,2)];                              const double C_1_0 = C_p[INDEX2(1,0,2)];
1838                              const register double C_0_1 = C_p[INDEX2(0,1,2)];                              const double C_0_1 = C_p[INDEX2(0,1,2)];
1839                              const register double C_1_1 = C_p[INDEX2(1,1,2)];                              const double C_1_1 = C_p[INDEX2(1,1,2)];
1840                              const register double C_0_2 = C_p[INDEX2(0,2,2)];                              const double C_0_2 = C_p[INDEX2(0,2,2)];
1841                              const register double C_1_2 = C_p[INDEX2(1,2,2)];                              const double C_1_2 = C_p[INDEX2(1,2,2)];
1842                              const register double C_0_3 = C_p[INDEX2(0,3,2)];                              const double C_0_3 = C_p[INDEX2(0,3,2)];
1843                              const register double C_1_3 = C_p[INDEX2(1,3,2)];                              const double C_1_3 = C_p[INDEX2(1,3,2)];
1844                              const register double tmp0_0 = C_1_0 + C_1_1;                              const double tmp0_0 = C_1_0 + C_1_1;
1845                              const register double tmp1_0 = C_1_2 + C_1_3;                              const double tmp1_0 = C_1_2 + C_1_3;
1846                              const register double tmp2_0 = C_0_1 + C_0_3;                              const double tmp2_0 = C_0_1 + C_0_3;
1847                              const register double tmp3_0 = C_0_0 + C_0_2;                              const double tmp3_0 = C_0_0 + C_0_2;
1848                              const register double tmp64_1 = C_0_2*w30;                              const double tmp64_1 = C_0_2*w30;
1849                              const register double tmp14_1 = C_0_2*w28;                              const double tmp14_1 = C_0_2*w28;
1850                              const register double tmp19_1 = C_0_3*w31;                              const double tmp19_1 = C_0_3*w31;
1851                              const register double tmp22_1 = C_1_0*w40;                              const double tmp22_1 = C_1_0*w40;
1852                              const register double tmp37_1 = tmp3_0*w35;                              const double tmp37_1 = tmp3_0*w35;
1853                              const register double tmp29_1 = C_0_1*w35;                              const double tmp29_1 = C_0_1*w35;
1854                              const register double tmp73_1 = C_0_2*w37;                              const double tmp73_1 = C_0_2*w37;
1855                              const register double tmp74_1 = C_1_2*w41;                              const double tmp74_1 = C_1_2*w41;
1856                              const register double tmp52_1 = C_1_3*w39;                              const double tmp52_1 = C_1_3*w39;
1857                              const register double tmp25_1 = C_1_1*w39;                              const double tmp25_1 = C_1_1*w39;
1858                              const register double tmp62_1 = C_0_1*w31;                              const double tmp62_1 = C_0_1*w31;
1859                              const register double tmp79_1 = C_1_1*w40;                              const double tmp79_1 = C_1_1*w40;
1860                              const register double tmp43_1 = C_1_1*w36;                              const double tmp43_1 = C_1_1*w36;
1861                              const register double tmp27_1 = C_0_3*w38;                              const double tmp27_1 = C_0_3*w38;
1862                              const register double tmp28_1 = C_1_0*w42;                              const double tmp28_1 = C_1_0*w42;
1863                              const register double tmp63_1 = C_1_1*w42;                              const double tmp63_1 = C_1_1*w42;
1864                              const register double tmp59_1 = C_0_3*w34;                              const double tmp59_1 = C_0_3*w34;
1865                              const register double tmp72_1 = C_1_3*w29;                              const double tmp72_1 = C_1_3*w29;
1866                              const register double tmp40_1 = tmp2_0*w35;                              const double tmp40_1 = tmp2_0*w35;
1867                              const register double tmp13_1 = C_0_3*w30;                              const double tmp13_1 = C_0_3*w30;
1868                              const register double tmp51_1 = C_1_2*w40;                              const double tmp51_1 = C_1_2*w40;
1869                              const register double tmp54_1 = C_1_2*w42;                              const double tmp54_1 = C_1_2*w42;
1870                              const register double tmp12_1 = C_0_0*w31;                              const double tmp12_1 = C_0_0*w31;
1871                              const register double tmp2_1 = tmp1_0*w32;                              const double tmp2_1 = tmp1_0*w32;
1872                              const register double tmp68_1 = C_0_2*w31;                              const double tmp68_1 = C_0_2*w31;
1873                              const register double tmp75_1 = C_1_0*w32;                              const double tmp75_1 = C_1_0*w32;
1874                              const register double tmp49_1 = C_1_1*w41;                              const double tmp49_1 = C_1_1*w41;
1875                              const register double tmp4_1 = C_0_2*w35;                              const double tmp4_1 = C_0_2*w35;
1876                              const register double tmp66_1 = C_0_3*w28;                              const double tmp66_1 = C_0_3*w28;
1877                              const register double tmp56_1 = C_0_1*w37;                              const double tmp56_1 = C_0_1*w37;
1878                              const register double tmp5_1 = C_0_1*w34;                              const double tmp5_1 = C_0_1*w34;
1879                              const register double tmp38_1 = tmp2_0*w34;                              const double tmp38_1 = tmp2_0*w34;
1880                              const register double tmp76_1 = C_0_1*w38;                              const double tmp76_1 = C_0_1*w38;
1881                              const register double tmp21_1 = C_0_1*w28;                              const double tmp21_1 = C_0_1*w28;
1882                              const register double tmp69_1 = C_0_1*w30;                              const double tmp69_1 = C_0_1*w30;
1883                              const register double tmp53_1 = C_1_0*w36;                              const double tmp53_1 = C_1_0*w36;
1884                              const register double tmp42_1 = C_1_2*w39;                              const double tmp42_1 = C_1_2*w39;
1885                              const register double tmp32_1 = tmp1_0*w29;                              const double tmp32_1 = tmp1_0*w29;
1886                              const register double tmp45_1 = C_1_0*w43;                              const double tmp45_1 = C_1_0*w43;
1887                              const register double tmp33_1 = tmp0_0*w32;                              const double tmp33_1 = tmp0_0*w32;
1888                              const register double tmp35_1 = C_1_0*w41;                              const double tmp35_1 = C_1_0*w41;
1889                              const register double tmp26_1 = C_1_2*w36;                              const double tmp26_1 = C_1_2*w36;
1890                              const register double tmp67_1 = C_0_0*w33;                              const double tmp67_1 = C_0_0*w33;
1891                              const register double tmp31_1 = C_0_2*w34;                              const double tmp31_1 = C_0_2*w34;
1892                              const register double tmp20_1 = C_0_0*w30;                              const double tmp20_1 = C_0_0*w30;
1893                              const register double tmp70_1 = C_0_0*w28;                              const double tmp70_1 = C_0_0*w28;
1894                              const register double tmp8_1 = tmp0_0*w39;                              const double tmp8_1 = tmp0_0*w39;
1895                              const register double tmp30_1 = C_1_3*w43;                              const double tmp30_1 = C_1_3*w43;
1896                              const register double tmp0_1 = tmp0_0*w29;                              const double tmp0_1 = tmp0_0*w29;
1897                              const register double tmp17_1 = C_1_3*w41;                              const double tmp17_1 = C_1_3*w41;
1898                              const register double tmp58_1 = C_0_0*w35;                              const double tmp58_1 = C_0_0*w35;
1899                              const register double tmp9_1 = tmp3_0*w33;                              const double tmp9_1 = tmp3_0*w33;
1900                              const register double tmp61_1 = C_1_3*w36;                              const double tmp61_1 = C_1_3*w36;
1901                              const register double tmp41_1 = tmp3_0*w34;                              const double tmp41_1 = tmp3_0*w34;
1902                              const register double tmp50_1 = C_1_3*w32;                              const double tmp50_1 = C_1_3*w32;
1903                              const register double tmp18_1 = C_1_1*w32;                              const double tmp18_1 = C_1_1*w32;
1904                              const register double tmp6_1 = tmp1_0*w36;                              const double tmp6_1 = tmp1_0*w36;
1905                              const register double tmp3_1 = C_0_0*w38;                              const double tmp3_1 = C_0_0*w38;
1906                              const register double tmp34_1 = C_1_1*w29;                              const double tmp34_1 = C_1_1*w29;
1907                              const register double tmp77_1 = C_0_3*w35;                              const double tmp77_1 = C_0_3*w35;
1908                              const register double tmp65_1 = C_1_2*w43;                              const double tmp65_1 = C_1_2*w43;
1909                              const register double tmp71_1 = C_0_3*w33;                              const double tmp71_1 = C_0_3*w33;
1910                              const register double tmp55_1 = C_1_1*w43;                              const double tmp55_1 = C_1_1*w43;
1911                              const register double tmp46_1 = tmp3_0*w28;                              const double tmp46_1 = tmp3_0*w28;
1912                              const register double tmp24_1 = C_0_0*w37;                              const double tmp24_1 = C_0_0*w37;
1913                              const register double tmp10_1 = tmp1_0*w39;                              const double tmp10_1 = tmp1_0*w39;
1914                              const register double tmp48_1 = C_1_0*w29;                              const double tmp48_1 = C_1_0*w29;
1915                              const register double tmp15_1 = C_0_1*w33;                              const double tmp15_1 = C_0_1*w33;
1916                              const register double tmp36_1 = C_1_2*w32;                              const double tmp36_1 = C_1_2*w32;
1917                              const register double tmp60_1 = C_1_0*w39;                              const double tmp60_1 = C_1_0*w39;
1918                              const register double tmp47_1 = tmp2_0*w33;                              const double tmp47_1 = tmp2_0*w33;
1919                              const register double tmp16_1 = C_1_2*w29;                              const double tmp16_1 = C_1_2*w29;
1920                              const register double tmp1_1 = C_0_3*w37;                              const double tmp1_1 = C_0_3*w37;
1921                              const register double tmp7_1 = tmp2_0*w28;                              const double tmp7_1 = tmp2_0*w28;
1922                              const register double tmp39_1 = C_1_3*w40;                              const double tmp39_1 = C_1_3*w40;
1923                              const register double tmp44_1 = C_1_3*w42;                              const double tmp44_1 = C_1_3*w42;
1924                              const register double tmp57_1 = C_0_2*w38;                              const double tmp57_1 = C_0_2*w38;
1925                              const register double tmp78_1 = C_0_0*w34;                              const double tmp78_1 = C_0_0*w34;
1926                              const register double tmp11_1 = tmp0_0*w36;                              const double tmp11_1 = tmp0_0*w36;
1927                              const register double tmp23_1 = C_0_2*w33;                              const double tmp23_1 = C_0_2*w33;
1928                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1929                              EM_S[INDEX2(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;
1930                              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 1942  void Rectangle::assemblePDESingle(Paso_S
1942                              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;
1943                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1944                          } else { // constant data                          } else { // constant data
1945                              const register double C_0 = C_p[0];                              const double C_0 = C_p[0];
1946                              const register double C_1 = C_p[1];                              const double C_1 = C_p[1];
1947                              const register double tmp0_1 = C_0*w47;                              const double tmp0_1 = C_0*w47;
1948                              const register double tmp1_1 = C_1*w45;                              const double tmp1_1 = C_1*w45;
1949                              const register double tmp2_1 = C_1*w48;                              const double tmp2_1 = C_1*w48;
1950                              const register double tmp3_1 = C_0*w51;                              const double tmp3_1 = C_0*w51;
1951                              const register double tmp4_1 = C_0*w44;                              const double tmp4_1 = C_0*w44;
1952                              const register double tmp5_1 = C_1*w49;                              const double tmp5_1 = C_1*w49;
1953                              const register double tmp6_1 = C_1*w50;                              const double tmp6_1 = C_1*w50;
1954                              const register double tmp7_1 = C_0*w46;                              const double tmp7_1 = C_0*w46;
1955                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1956                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
1957                              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 1977  void Rectangle::assemblePDESingle(Paso_S
1977                          add_EM_S=true;                          add_EM_S=true;
1978                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
1979                          if (D.actsExpanded()) {                          if (D.actsExpanded()) {
1980                              const register double D_0 = D_p[0];                              const double D_0 = D_p[0];
1981                              const register double D_1 = D_p[1];                              const double D_1 = D_p[1];
1982                              const register double D_2 = D_p[2];                              const double D_2 = D_p[2];
1983                              const register double D_3 = D_p[3];                              const double D_3 = D_p[3];
1984                              const register double tmp0_0 = D_0 + D_1;                              const double tmp0_0 = D_0 + D_1;
1985                              const register double tmp1_0 = D_2 + D_3;                              const double tmp1_0 = D_2 + D_3;
1986                              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;
1987                              const register double tmp3_0 = D_1 + D_2;                              const double tmp3_0 = D_1 + D_2;
1988                              const register double tmp4_0 = D_1 + D_3;                              const double tmp4_0 = D_1 + D_3;
1989                              const register double tmp5_0 = D_0 + D_2;                              const double tmp5_0 = D_0 + D_2;
1990                              const register double tmp6_0 = D_0 + D_3;                              const double tmp6_0 = D_0 + D_3;
1991                              const register double tmp0_1 = tmp0_0*w52;                              const double tmp0_1 = tmp0_0*w52;
1992                              const register double tmp1_1 = tmp1_0*w53;                              const double tmp1_1 = tmp1_0*w53;
1993                              const register double tmp2_1 = tmp2_0*w54;                              const double tmp2_1 = tmp2_0*w54;
1994                              const register double tmp3_1 = tmp1_0*w52;                              const double tmp3_1 = tmp1_0*w52;
1995                              const register double tmp4_1 = tmp0_0*w53;                              const double tmp4_1 = tmp0_0*w53;
1996                              const register double tmp5_1 = tmp3_0*w54;                              const double tmp5_1 = tmp3_0*w54;
1997                              const register double tmp6_1 = D_0*w55;                              const double tmp6_1 = D_0*w55;
1998                              const register double tmp7_1 = D_3*w56;                              const double tmp7_1 = D_3*w56;
1999                              const register double tmp8_1 = D_3*w55;                              const double tmp8_1 = D_3*w55;
2000                              const register double tmp9_1 = D_0*w56;                              const double tmp9_1 = D_0*w56;
2001                              const register double tmp10_1 = tmp4_0*w52;                              const double tmp10_1 = tmp4_0*w52;
2002                              const register double tmp11_1 = tmp5_0*w53;                              const double tmp11_1 = tmp5_0*w53;
2003                              const register double tmp12_1 = tmp5_0*w52;                              const double tmp12_1 = tmp5_0*w52;
2004                              const register double tmp13_1 = tmp4_0*w53;                              const double tmp13_1 = tmp4_0*w53;
2005                              const register double tmp14_1 = tmp6_0*w54;                              const double tmp14_1 = tmp6_0*w54;
2006                              const register double tmp15_1 = D_2*w55;                              const double tmp15_1 = D_2*w55;
2007                              const register double tmp16_1 = D_1*w56;                              const double tmp16_1 = D_1*w56;
2008                              const register double tmp17_1 = D_1*w55;                              const double tmp17_1 = D_1*w55;
2009                              const register double tmp18_1 = D_2*w56;                              const double tmp18_1 = D_2*w56;
2010                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2011                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2012                              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 2024  void Rectangle::assemblePDESingle(Paso_S
2024                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2025                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2026                          } else { // constant data                          } else { // constant data
2027                              const register double tmp0_1 = D_p[0]*w57;                              const double tmp0_1 = D_p[0]*w57;
2028                              const register double tmp1_1 = D_p[0]*w58;                              const double tmp1_1 = D_p[0]*w58;
2029                              const register double tmp2_1 = D_p[0]*w59;                              const double tmp2_1 = D_p[0]*w59;
2030                              EM_S[INDEX2(0,0,4)]+=tmp2_1;                              EM_S[INDEX2(0,0,4)]+=tmp2_1;
2031                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
2032                              EM_S[INDEX2(2,0,4)]+=tmp0_1;                              EM_S[INDEX2(2,0,4)]+=tmp0_1;
# Line 2066  void Rectangle::assemblePDESingle(Paso_S Line 2052  void Rectangle::assemblePDESingle(Paso_S
2052                          add_EM_F=true;                          add_EM_F=true;
2053                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2054                          if (X.actsExpanded()) {                          if (X.actsExpanded()) {
2055                              const register double X_0_0 = X_p[INDEX2(0,0,2)];                              const double X_0_0 = X_p[INDEX2(0,0,2)];
2056                              const register double X_1_0 = X_p[INDEX2(1,0,2)];                              const double X_1_0 = X_p[INDEX2(1,0,2)];
2057                              const register double X_0_1 = X_p[INDEX2(0,1,2)];                              const double X_0_1 = X_p[INDEX2(0,1,2)];
2058                              const register double X_1_1 = X_p[INDEX2(1,1,2)];                              const double X_1_1 = X_p[INDEX2(1,1,2)];
2059                              const register double X_0_2 = X_p[INDEX2(0,2,2)];                              const double X_0_2 = X_p[INDEX2(0,2,2)];
2060                              const register double X_1_2 = X_p[INDEX2(1,2,2)];                              const double X_1_2 = X_p[INDEX2(1,2,2)];
2061                              const register double X_0_3 = X_p[INDEX2(0,3,2)];                              const double X_0_3 = X_p[INDEX2(0,3,2)];
2062                              const register double X_1_3 = X_p[INDEX2(1,3,2)];                              const double X_1_3 = X_p[INDEX2(1,3,2)];
2063                              const register double tmp0_0 = X_0_2 + X_0_3;                              const double tmp0_0 = X_0_2 + X_0_3;
2064                              const register double tmp1_0 = X_1_1 + X_1_3;                              const double tmp1_0 = X_1_1 + X_1_3;
2065                              const register double tmp2_0 = X_1_0 + X_1_2;                              const double tmp2_0 = X_1_0 + X_1_2;
2066                              const register double tmp3_0 = X_0_0 + X_0_1;                              const double tmp3_0 = X_0_0 + X_0_1;
2067                              const register double tmp0_1 = tmp0_0*w63;                              const double tmp0_1 = tmp0_0*w63;
2068                              const register double tmp1_1 = tmp1_0*w62;                              const double tmp1_1 = tmp1_0*w62;
2069                              const register double tmp2_1 = tmp2_0*w61;                              const double tmp2_1 = tmp2_0*w61;
2070                              const register double tmp3_1 = tmp3_0*w60;                              const double tmp3_1 = tmp3_0*w60;
2071                              const register double tmp4_1 = tmp0_0*w65;                              const double tmp4_1 = tmp0_0*w65;
2072                              const register double tmp5_1 = tmp3_0*w64;                              const double tmp5_1 = tmp3_0*w64;
2073                              const register double tmp6_1 = tmp2_0*w62;                              const double tmp6_1 = tmp2_0*w62;
2074                              const register double tmp7_1 = tmp1_0*w61;                              const double tmp7_1 = tmp1_0*w61;
2075                              const register double tmp8_1 = tmp2_0*w66;                              const double tmp8_1 = tmp2_0*w66;
2076                              const register double tmp9_1 = tmp3_0*w63;                              const double tmp9_1 = tmp3_0*w63;
2077                              const register double tmp10_1 = tmp0_0*w60;                              const double tmp10_1 = tmp0_0*w60;
2078                              const register double tmp11_1 = tmp1_0*w67;                              const double tmp11_1 = tmp1_0*w67;
2079                              const register double tmp12_1 = tmp1_0*w66;                              const double tmp12_1 = tmp1_0*w66;
2080                              const register double tmp13_1 = tmp3_0*w65;                              const double tmp13_1 = tmp3_0*w65;
2081                              const register double tmp14_1 = tmp0_0*w64;                              const double tmp14_1 = tmp0_0*w64;
2082                              const register double tmp15_1 = tmp2_0*w67;                              const double tmp15_1 = tmp2_0*w67;
2083                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2084                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2085                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2086                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2087                          } else { // constant data                          } else { // constant data
2088                              const register double tmp0_1 = X_p[1]*w69;                              const double tmp0_1 = X_p[1]*w69;
2089                              const register double tmp1_1 = X_p[0]*w68;                              const double tmp1_1 = X_p[0]*w68;
2090                              const register double tmp2_1 = X_p[0]*w70;                              const double tmp2_1 = X_p[0]*w70;
2091                              const register double tmp3_1 = X_p[1]*w71;                              const double tmp3_1 = X_p[1]*w71;
2092                              EM_F[0]+=tmp0_1 + tmp1_1;                              EM_F[0]+=tmp0_1 + tmp1_1;
2093                              EM_F[1]+=tmp0_1 + tmp2_1;                              EM_F[1]+=tmp0_1 + tmp2_1;
2094                              EM_F[2]+=tmp1_1 + tmp3_1;                              EM_F[2]+=tmp1_1 + tmp3_1;
# Line 2116  void Rectangle::assemblePDESingle(Paso_S Line 2102  void Rectangle::assemblePDESingle(Paso_S
2102                          add_EM_F=true;                          add_EM_F=true;
2103                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2104                          if (Y.actsExpanded()) {                          if (Y.actsExpanded()) {
2105                              const register double Y_0 = Y_p[0];                              const double Y_0 = Y_p[0];
2106                              const register double Y_1 = Y_p[1];                              const double Y_1 = Y_p[1];
2107                              const register double Y_2 = Y_p[2];                              const double Y_2 = Y_p[2];
2108                              const register double Y_3 = Y_p[3];                              const double Y_3 = Y_p[3];
2109                              const register double tmp0_0 = Y_1 + Y_2;                              const double tmp0_0 = Y_1 + Y_2;
2110                              const register double tmp1_0 = Y_0 + Y_3;                              const double tmp1_0 = Y_0 + Y_3;
2111                              const register double tmp0_1 = Y_0*w72;                              const double tmp0_1 = Y_0*w72;
2112                              const register double tmp1_1 = tmp0_0*w73;                              const double tmp1_1 = tmp0_0*w73;
2113                              const register double tmp2_1 = Y_3*w74;                              const double tmp2_1 = Y_3*w74;
2114                              const register double tmp3_1 = Y_1*w72;                              const double tmp3_1 = Y_1*w72;
2115                              const register double tmp4_1 = tmp1_0*w73;                              const double tmp4_1 = tmp1_0*w73;
2116                              const register double tmp5_1 = Y_2*w74;                              const double tmp5_1 = Y_2*w74;
2117                              const register double tmp6_1 = Y_2*w72;                              const double tmp6_1 = Y_2*w72;
2118                              const register double tmp7_1 = Y_1*w74;                              const double tmp7_1 = Y_1*w74;
2119                              const register double tmp8_1 = Y_3*w72;                              const double tmp8_1 = Y_3*w72;
2120                              const register double tmp9_1 = Y_0*w74;                              const double tmp9_1 = Y_0*w74;
2121                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2122                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2123                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2124                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2125                          } else { // constant data                          } else { // constant data
2126                              const register double tmp0_1 = Y_p[0]*w75;                              const double tmp0_1 = Y_p[0]*w75;
2127                              EM_F[0]+=tmp0_1;                              EM_F[0]+=tmp0_1;
2128                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
2129                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
# Line 2174  void Rectangle::assemblePDESingle(Paso_S Line 2160  void Rectangle::assemblePDESingle(Paso_S
2160  void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2161          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2162          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2163          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  
2164  {  {
2165      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
2166      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 2214  void Rectangle::assemblePDESingleReduced Line 2199  void Rectangle::assemblePDESingleReduced
2199                      if (!A.isEmpty()) {                      if (!A.isEmpty()) {
2200                          add_EM_S=true;                          add_EM_S=true;
2201                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2202                          const register double A_00 = A_p[INDEX2(0,0,2)];                          const double A_00 = A_p[INDEX2(0,0,2)];
2203                          const register double A_10 = A_p[INDEX2(1,0,2)];                          const double A_10 = A_p[INDEX2(1,0,2)];
2204                          const register double A_01 = A_p[INDEX2(0,1,2)];                          const double A_01 = A_p[INDEX2(0,1,2)];
2205                          const register double A_11 = A_p[INDEX2(1,1,2)];                          const double A_11 = A_p[INDEX2(1,1,2)];
2206                          const register double tmp0_0 = A_01 + A_10;                          const double tmp0_0 = A_01 + A_10;
2207                          const register double tmp0_1 = A_11*w3;                          const double tmp0_1 = A_11*w3;
2208                          const register double tmp1_1 = A_00*w0;                          const double tmp1_1 = A_00*w0;
2209                          const register double tmp2_1 = A_01*w1;                          const double tmp2_1 = A_01*w1;
2210                          const register double tmp3_1 = A_10*w2;                          const double tmp3_1 = A_10*w2;
2211                          const register double tmp4_1 = tmp0_0*w1;                          const double tmp4_1 = tmp0_0*w1;
2212                          const register double tmp5_1 = A_11*w4;                          const double tmp5_1 = A_11*w4;
2213                          const register double tmp6_1 = A_00*w5;                          const double tmp6_1 = A_00*w5;
2214                          const register double tmp7_1 = tmp0_0*w2;                          const double tmp7_1 = tmp0_0*w2;
2215                          const register double tmp8_1 = A_10*w1;                          const double tmp8_1 = A_10*w1;
2216                          const register double tmp9_1 = A_01*w2;                          const double tmp9_1 = A_01*w2;
2217                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2218                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2219                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
# Line 2252  void Rectangle::assemblePDESingleReduced Line 2237  void Rectangle::assemblePDESingleReduced
2237                      if (!B.isEmpty()) {                      if (!B.isEmpty()) {
2238                          add_EM_S=true;                          add_EM_S=true;
2239                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2240                          const register double tmp2_1 = B_p[0]*w8;                          const double tmp2_1 = B_p[0]*w8;
2241                          const register double tmp0_1 = B_p[0]*w6;                          const double tmp0_1 = B_p[0]*w6;
2242                          const register double tmp3_1 = B_p[1]*w9;                          const double tmp3_1 = B_p[1]*w9;
2243                          const register double tmp1_1 = B_p[1]*w7;                          const double tmp1_1 = B_p[1]*w7;
2244                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2245                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2246                          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 2264  void Rectangle::assemblePDESingleReduced
2264                      if (!C.isEmpty()) {                      if (!C.isEmpty()) {
2265                          add_EM_S=true;                          add_EM_S=true;
2266                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2267                          const register double tmp1_1 = C_p[1]*w7;                          const double tmp1_1 = C_p[1]*w7;
2268                          const register double tmp0_1 = C_p[0]*w8;                          const double tmp0_1 = C_p[0]*w8;
2269                          const register double tmp3_1 = C_p[0]*w6;                          const double tmp3_1 = C_p[0]*w6;
2270                          const register double tmp2_1 = C_p[1]*w9;                          const double tmp2_1 = C_p[1]*w9;
2271                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2272                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2273                          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 2291  void Rectangle::assemblePDESingleReduced
2291                      if (!D.isEmpty()) {                      if (!D.isEmpty()) {
2292                          add_EM_S=true;                          add_EM_S=true;
2293                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2294                          const register double tmp0_1 = D_p[0]*w10;                          const double tmp0_1 = D_p[0]*w10;
2295                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1;
2296                          EM_S[INDEX2(1,0,4)]+=tmp0_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1;
2297                          EM_S[INDEX2(2,0,4)]+=tmp0_1;                          EM_S[INDEX2(2,0,4)]+=tmp0_1;
# Line 2330  void Rectangle::assemblePDESingleReduced Line 2315  void Rectangle::assemblePDESingleReduced
2315                      if (!X.isEmpty()) {                      if (!X.isEmpty()) {
2316                          add_EM_F=true;                          add_EM_F=true;
2317                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2318                          const register double tmp0_1 = X_p[0]*w11;                          const double tmp0_1 = X_p[0]*w11;
2319                          const register double tmp2_1 = X_p[0]*w13;                          const double tmp2_1 = X_p[0]*w13;
2320                          const register double tmp1_1 = X_p[1]*w12;                          const double tmp1_1 = X_p[1]*w12;
2321                          const register double tmp3_1 = X_p[1]*w14;                          const double tmp3_1 = X_p[1]*w14;
2322                          EM_F[0]+=tmp0_1 + tmp1_1;                          EM_F[0]+=tmp0_1 + tmp1_1;
2323                          EM_F[1]+=tmp1_1 + tmp2_1;                          EM_F[1]+=tmp1_1 + tmp2_1;
2324                          EM_F[2]+=tmp0_1 + tmp3_1;                          EM_F[2]+=tmp0_1 + tmp3_1;
# Line 2345  void Rectangle::assemblePDESingleReduced Line 2330  void Rectangle::assemblePDESingleReduced
2330                      if (!Y.isEmpty()) {                      if (!Y.isEmpty()) {
2331                          add_EM_F=true;                          add_EM_F=true;
2332                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2333                          const register double tmp0_1 = Y_p[0]*w15;                          const double tmp0_1 = Y_p[0]*w15;
2334                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0_1;
2335                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0_1;
2336                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
# Line 2381  void Rectangle::assemblePDESingleReduced Line 2366  void Rectangle::assemblePDESingleReduced
2366  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2367          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2368          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2369          const escript::Data& X, const escript::Data& Y,          const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
2370  {  {
2371      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
2372      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 2492  void Rectangle::assemblePDESystem(Paso_S Line 2476  void Rectangle::assemblePDESystem(Paso_S
2476                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
2477                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2478                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2479                                      const register double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];                                      const double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];
2480                                      const register double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];                                      const double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];
2481                                      const register double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];                                      const double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];
2482                                      const register double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];                                      const double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];
2483                                      const register double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];                                      const double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];
2484                                      const register double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];                                      const double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];
2485                                      const register double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];                                      const double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];
2486                                      const register double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];                                      const double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];
2487                                      const register double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];                                      const double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];
2488                                      const register double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];                                      const double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];
2489                                      const register double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];                                      const double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];
2490                                      const register double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];                                      const double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];
2491                                      const register double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];                                      const double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];
2492                                      const register double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];                                      const double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];
2493                                      const register double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];                                      const double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];
2494                                      const register double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];                                      const double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];
2495                                      const register double tmp0_0 = A_01_0 + A_01_3;                                      const double tmp0_0 = A_01_0 + A_01_3;
2496                                      const register double tmp1_0 = A_00_0 + A_00_1;                                      const double tmp1_0 = A_00_0 + A_00_1;
2497                                      const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                                      const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
2498                                      const register double tmp3_0 = A_00_2 + A_00_3;                                      const double tmp3_0 = A_00_2 + A_00_3;
2499                                      const register double tmp4_0 = A_10_1 + A_10_2;                                      const double tmp4_0 = A_10_1 + A_10_2;
2500                                      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;
2501                                      const register double tmp6_0 = A_01_3 + A_10_0;                                      const double tmp6_0 = A_01_3 + A_10_0;
2502                                      const register double tmp7_0 = A_01_0 + A_10_3;                                      const double tmp7_0 = A_01_0 + A_10_3;
2503                                      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;
2504                                      const register double tmp9_0 = A_01_0 + A_10_0;                                      const double tmp9_0 = A_01_0 + A_10_0;
2505                                      const register double tmp10_0 = A_01_3 + A_10_3;                                      const double tmp10_0 = A_01_3 + A_10_3;
2506                                      const register double tmp11_0 = A_11_1 + A_11_3;                                      const double tmp11_0 = A_11_1 + A_11_3;
2507                                      const register double tmp12_0 = A_11_0 + A_11_2;                                      const double tmp12_0 = A_11_0 + A_11_2;
2508                                      const register double tmp13_0 = A_01_2 + A_10_1;                                      const double tmp13_0 = A_01_2 + A_10_1;
2509                                      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;
2510                                      const register double tmp15_0 = A_01_1 + A_10_2;                                      const double tmp15_0 = A_01_1 + A_10_2;
2511                                      const register double tmp16_0 = A_10_0 + A_10_3;                                      const double tmp16_0 = A_10_0 + A_10_3;
2512                                      const register double tmp17_0 = A_01_1 + A_01_2;                                      const double tmp17_0 = A_01_1 + A_01_2;
2513                                      const register double tmp18_0 = A_01_1 + A_10_1;                                      const double tmp18_0 = A_01_1 + A_10_1;
2514                                      const register double tmp19_0 = A_01_2 + A_10_2;                                      const double tmp19_0 = A_01_2 + A_10_2;
2515                                      const register double tmp0_1 = A_10_3*w8;                                      const double tmp0_1 = A_10_3*w8;
2516                                      const register double tmp1_1 = tmp0_0*w1;                                      const double tmp1_1 = tmp0_0*w1;
2517                                      const register double tmp2_1 = A_01_1*w4;                                      const double tmp2_1 = A_01_1*w4;
2518                                      const register double tmp3_1 = tmp1_0*w0;                                      const double tmp3_1 = tmp1_0*w0;
2519                                      const register double tmp4_1 = A_01_2*w7;                                      const double tmp4_1 = A_01_2*w7;
2520                                      const register double tmp5_1 = tmp2_0*w3;                                      const double tmp5_1 = tmp2_0*w3;
2521                                      const register double tmp6_1 = tmp3_0*w6;                                      const double tmp6_1 = tmp3_0*w6;
2522                                      const register double tmp7_1 = A_10_0*w2;                                      const double tmp7_1 = A_10_0*w2;
2523                                      const register double tmp8_1 = tmp4_0*w5;                                      const double tmp8_1 = tmp4_0*w5;
2524                                      const register double tmp9_1 = tmp2_0*w10;                                      const double tmp9_1 = tmp2_0*w10;
2525                                      const register double tmp10_1 = tmp5_0*w9;                                      const double tmp10_1 = tmp5_0*w9;
2526                                      const register double tmp11_1 = tmp6_0*w7;                                      const double tmp11_1 = tmp6_0*w7;
2527                                      const register double tmp12_1 = tmp7_0*w4;                                      const double tmp12_1 = tmp7_0*w4;
2528                                      const register double tmp13_1 = tmp8_0*w1;                                      const double tmp13_1 = tmp8_0*w1;
2529                                      const register double tmp14_1 = A_10_0*w8;                                      const double tmp14_1 = A_10_0*w8;
2530                                      const register double tmp15_1 = A_01_2*w4;                                      const double tmp15_1 = A_01_2*w4;
2531                                      const register double tmp16_1 = tmp3_0*w0;                                      const double tmp16_1 = tmp3_0*w0;
2532                                      const register double tmp17_1 = A_01_1*w7;                                      const double tmp17_1 = A_01_1*w7;
2533                                      const register double tmp18_1 = tmp1_0*w6;                                      const double tmp18_1 = tmp1_0*w6;
2534                                      const register double tmp19_1 = A_10_3*w2;                                      const double tmp19_1 = A_10_3*w2;
2535                                      const register double tmp20_1 = tmp9_0*w4;                                      const double tmp20_1 = tmp9_0*w4;
2536                                      const register double tmp21_1 = tmp1_0*w11;                                      const double tmp21_1 = tmp1_0*w11;
2537                                      const register double tmp22_1 = tmp10_0*w7;                                      const double tmp22_1 = tmp10_0*w7;
2538                                      const register double tmp23_1 = tmp3_0*w14;                                      const double tmp23_1 = tmp3_0*w14;
2539                                      const register double tmp24_1 = tmp11_0*w13;                                      const double tmp24_1 = tmp11_0*w13;
2540                                      const register double tmp25_1 = tmp12_0*w12;                                      const double tmp25_1 = tmp12_0*w12;
2541                                      const register double tmp26_1 = tmp10_0*w4;                                      const double tmp26_1 = tmp10_0*w4;
2542                                      const register double tmp27_1 = tmp3_0*w11;                                      const double tmp27_1 = tmp3_0*w11;
2543                                      const register double tmp28_1 = tmp9_0*w7;                                      const double tmp28_1 = tmp9_0*w7;
2544                                      const register double tmp29_1 = tmp1_0*w14;                                      const double tmp29_1 = tmp1_0*w14;
2545                                      const register double tmp30_1 = tmp12_0*w13;                                      const double tmp30_1 = tmp12_0*w13;
2546                                      const register double tmp31_1 = tmp11_0*w12;                                      const double tmp31_1 = tmp11_0*w12;
2547                                      const register double tmp32_1 = tmp13_0*w2;                                      const double tmp32_1 = tmp13_0*w2;
2548                                      const register double tmp33_1 = tmp14_0*w5;                                      const double tmp33_1 = tmp14_0*w5;
2549                                      const register double tmp34_1 = tmp15_0*w8;                                      const double tmp34_1 = tmp15_0*w8;
2550                                      const register double tmp35_1 = A_01_0*w8;                                      const double tmp35_1 = A_01_0*w8;
2551                                      const register double tmp36_1 = tmp16_0*w1;                                      const double tmp36_1 = tmp16_0*w1;
2552                                      const register double tmp37_1 = A_10_1*w4;                                      const double tmp37_1 = A_10_1*w4;
2553                                      const register double tmp38_1 = tmp5_0*w15;                                      const double tmp38_1 = tmp5_0*w15;
2554                                      const register double tmp39_1 = A_10_2*w7;                                      const double tmp39_1 = A_10_2*w7;
2555                                      const register double tmp40_1 = tmp11_0*w17;                                      const double tmp40_1 = tmp11_0*w17;
2556                                      const register double tmp41_1 = A_01_3*w2;                                      const double tmp41_1 = A_01_3*w2;
2557                                      const register double tmp42_1 = tmp12_0*w16;                                      const double tmp42_1 = tmp12_0*w16;
2558                                      const register double tmp43_1 = tmp17_0*w5;                                      const double tmp43_1 = tmp17_0*w5;
2559                                      const register double tmp44_1 = tmp7_0*w7;                                      const double tmp44_1 = tmp7_0*w7;
2560                                      const register double tmp45_1 = tmp6_0*w4;                                      const double tmp45_1 = tmp6_0*w4;
2561                                      const register double tmp46_1 = A_01_3*w8;                                      const double tmp46_1 = A_01_3*w8;
2562                                      const register double tmp47_1 = A_10_2*w4;                                      const double tmp47_1 = A_10_2*w4;
2563                                      const register double tmp48_1 = A_10_1*w7;                                      const double tmp48_1 = A_10_1*w7;
2564                                      const register double tmp49_1 = tmp12_0*w17;                                      const double tmp49_1 = tmp12_0*w17;
2565                                      const register double tmp50_1 = A_01_0*w2;                                      const double tmp50_1 = A_01_0*w2;
2566                                      const register double tmp51_1 = tmp11_0*w16;                                      const double tmp51_1 = tmp11_0*w16;
2567                                      const register double tmp52_1 = tmp18_0*w8;                                      const double tmp52_1 = tmp18_0*w8;
2568                                      const register double tmp53_1 = tmp19_0*w2;                                      const double tmp53_1 = tmp19_0*w2;
2569                                      const register double tmp54_1 = tmp13_0*w8;                                      const double tmp54_1 = tmp13_0*w8;
2570                                      const register double tmp55_1 = tmp15_0*w2;                                      const double tmp55_1 = tmp15_0*w2;
2571                                      const register double tmp56_1 = tmp19_0*w8;                                      const double tmp56_1 = tmp19_0*w8;
2572                                      const register double tmp57_1 = tmp18_0*w2;                                      const double tmp57_1 = tmp18_0*w2;
2573                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
2574                                      EM_S[INDEX4(k,m,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;
2575                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;                                      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 2591  void Rectangle::assemblePDESystem(Paso_S
2591                          } else { // constant data                          } else { // constant data
2592                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2593                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2594                                      const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];                                      const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
2595                                      const register double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];                                      const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
2596                                      const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];                                      const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
2597                                      const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];                                      const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
2598                                      const register double tmp0_0 = A_01 + A_10;                                      const double tmp0_0 = A_01 + A_10;
2599                                      const register double tmp0_1 = A_00*w18;                                      const double tmp0_1 = A_00*w18;
2600                                      const register double tmp1_1 = A_01*w19;                                      const double tmp1_1 = A_01*w19;
2601                                      const register double tmp2_1 = A_10*w20;                                      const double tmp2_1 = A_10*w20;
2602                                      const register double tmp3_1 = A_11*w21;                                      const double tmp3_1 = A_11*w21;
2603                                      const register double tmp4_1 = A_00*w22;                                      const double tmp4_1 = A_00*w22;
2604                                      const register double tmp5_1 = tmp0_0*w19;                                      const double tmp5_1 = tmp0_0*w19;
2605                                      const register double tmp6_1 = A_11*w23;                                      const double tmp6_1 = A_11*w23;
2606                                      const register double tmp7_1 = A_11*w25;                                      const double tmp7_1 = A_11*w25;
2607                                      const register double tmp8_1 = A_00*w24;                                      const double tmp8_1 = A_00*w24;
2608                                      const register double tmp9_1 = tmp0_0*w20;                                      const double tmp9_1 = tmp0_0*w20;
2609                                      const register double tmp10_1 = A_01*w20;                                      const double tmp10_1 = A_01*w20;
2610                                      const register double tmp11_1 = A_11*w27;                                      const double tmp11_1 = A_11*w27;
2611                                      const register double tmp12_1 = A_00*w26;                                      const double tmp12_1 = A_00*w26;
2612                                      const register double tmp13_1 = A_10*w19;                                      const double tmp13_1 = A_10*w19;
2613                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2614                                      EM_S[INDEX4(k,m,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;
2615                                      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 2639  void Rectangle::assemblePDESystem(Paso_S
2639                          if (B.actsExpanded()) {                          if (B.actsExpanded()) {
2640                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2641                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2642                                      const register double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];                                      const double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
2643                                      const register double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];                                      const double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
2644                                      const register double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];                                      const double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
2645                                      const register double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];                                      const double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
2646                                      const register double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];                                      const double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];
2647                                      const register double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];                                      const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
2648                                      const register double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];                                      const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
2649                                      const register double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];                                      const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
2650                                      const register double tmp0_0 = B_1_0 + B_1_1;                                      const double tmp0_0 = B_1_0 + B_1_1;
2651                                      const register double tmp1_0 = B_1_2 + B_1_3;                                      const double tmp1_0 = B_1_2 + B_1_3;
2652                                      const register double tmp2_0 = B_0_1 + B_0_3;                                      const double tmp2_0 = B_0_1 + B_0_3;
2653                                      const register double tmp3_0 = B_0_0 + B_0_2;                                      const double tmp3_0 = B_0_0 + B_0_2;
2654                                      const register double tmp63_1 = B_1_1*w42;                                      const double tmp63_1 = B_1_1*w42;
2655                                      const register double tmp79_1 = B_1_1*w40;                                      const double tmp79_1 = B_1_1*w40;
2656                                      const register double tmp37_1 = tmp3_0*w35;                                      const double tmp37_1 = tmp3_0*w35;
2657                                      const register double tmp8_1 = tmp0_0*w32;                                      const double tmp8_1 = tmp0_0*w32;
2658                                      const register double tmp71_1 = B_0_1*w34;                                      const double tmp71_1 = B_0_1*w34;
2659                                      const register double tmp19_1 = B_0_3*w31;                                      const double tmp19_1 = B_0_3*w31;
2660                                      const register double tmp15_1 = B_0_3*w34;                                      const double tmp15_1 = B_0_3*w34;
2661                                      const register double tmp9_1 = tmp3_0*w34;                                      const double tmp9_1 = tmp3_0*w34;
2662                                      const register double tmp35_1 = B_1_0*w36;                                      const double tmp35_1 = B_1_0*w36;
2663                                      const register double tmp66_1 = B_0_3*w28;                                      const double tmp66_1 = B_0_3*w28;
2664                                      const register double tmp28_1 = B_1_0*w42;                                      const double tmp28_1 = B_1_0*w42;
2665                                      const register double tmp22_1 = B_1_0*w40;                                      const double tmp22_1 = B_1_0*w40;
2666                                      const register double tmp16_1 = B_1_2*w29;                                      const double tmp16_1 = B_1_2*w29;
2667                                      const register double tmp6_1 = tmp2_0*w35;                                      const double tmp6_1 = tmp2_0*w35;
2668                                      const register double tmp55_1 = B_1_3*w40;                                      const double tmp55_1 = B_1_3*w40;
2669                                      const register double tmp50_1 = B_1_3*w42;                                      const double tmp50_1 = B_1_3*w42;
2670                                      const register double tmp7_1 = tmp1_0*w29;                                      const double tmp7_1 = tmp1_0*w29;
2671                                      const register double tmp1_1 = tmp1_0*w32;                                      const double tmp1_1 = tmp1_0*w32;
2672                                      const register double tmp57_1 = B_0_3*w30;                                      const double tmp57_1 = B_0_3*w30;
2673                                      const register double tmp18_1 = B_1_1*w32;                                      const double tmp18_1 = B_1_1*w32;
2674                                      const register double tmp53_1 = B_1_0*w41;                                      const double tmp53_1 = B_1_0*w41;
2675                                      const register double tmp61_1 = B_1_3*w36;                                      const double tmp61_1 = B_1_3*w36;
2676                                      const register double tmp27_1 = B_0_3*w38;                                      const double tmp27_1 = B_0_3*w38;
2677                                      const register double tmp64_1 = B_0_2*w30;                                      const double tmp64_1 = B_0_2*w30;
2678                                      const register double tmp76_1 = B_0_1*w38;                                      const double tmp76_1 = B_0_1*w38;
2679                                      const register double tmp39_1 = tmp2_0*w34;                                      const double tmp39_1 = tmp2_0*w34;
2680                                      const register double tmp62_1 = B_0_1*w31;                                      const double tmp62_1 = B_0_1*w31;
2681                                      const register double tmp56_1 = B_0_0*w31;                                      const double tmp56_1 = B_0_0*w31;
2682                                      const register double tmp49_1 = B_1_1*w36;                                      const double tmp49_1 = B_1_1*w36;
2683                                      const register double tmp2_1 = B_0_2*w31;                                      const double tmp2_1 = B_0_2*w31;
2684                                      const register double tmp23_1 = B_0_2*w33;                                      const double tmp23_1 = B_0_2*w33;
2685                                      const register double tmp38_1 = B_1_1*w43;                                      const double tmp38_1 = B_1_1*w43;
2686                                      const register double tmp74_1 = B_1_2*w41;                                      const double tmp74_1 = B_1_2*w41;
2687                                      const register double tmp43_1 = B_1_1*w41;                                      const double tmp43_1 = B_1_1*w41;
2688                                      const register double tmp58_1 = B_0_2*w28;                                      const double tmp58_1 = B_0_2*w28;
2689                                      const register double tmp67_1 = B_0_0*w33;                                      const double tmp67_1 = B_0_0*w33;
2690                                      const register double tmp33_1 = tmp0_0*w39;                                      const double tmp33_1 = tmp0_0*w39;
2691                                      const register double tmp4_1 = B_0_0*w28;                                      const double tmp4_1 = B_0_0*w28;
2692                                      const register double tmp20_1 = B_0_0*w30;                                      const double tmp20_1 = B_0_0*w30;
2693                                      const register double tmp13_1 = B_0_2*w38;                                      const double tmp13_1 = B_0_2*w38;
2694                                      const register double tmp65_1 = B_1_2*w43;                                      const double tmp65_1 = B_1_2*w43;
2695                                      const register double tmp0_1 = tmp0_0*w29;                                      const double tmp0_1 = tmp0_0*w29;
2696                                      const register double tmp41_1 = tmp3_0*w33;                                      const double tmp41_1 = tmp3_0*w33;
2697                                      const register double tmp73_1 = B_0_2*w37;                                      const double tmp73_1 = B_0_2*w37;
2698                                      const register double tmp69_1 = B_0_0*w38;                                      const double tmp69_1 = B_0_0*w38;
2699                                      const register double tmp48_1 = B_1_2*w39;                                      const double tmp48_1 = B_1_2*w39;
2700                                      const register double tmp59_1 = B_0_1*w33;                                      const double tmp59_1 = B_0_1*w33;
2701                                      const register double tmp17_1 = B_1_3*w41;                                      const double tmp17_1 = B_1_3*w41;
2702                                      const register double tmp5_1 = B_0_3*w33;                                      const double tmp5_1 = B_0_3*w33;
2703                                      const register double tmp3_1 = B_0_1*w30;                                      const double tmp3_1 = B_0_1*w30;
2704                                      const register double tmp21_1 = B_0_1*w28;                                      const double tmp21_1 = B_0_1*w28;
2705                                      const register double tmp42_1 = B_1_0*w29;                                      const double tmp42_1 = B_1_0*w29;
2706                                      const register double tmp54_1 = B_1_2*w32;                                      const double tmp54_1 = B_1_2*w32;
2707                                      const register double tmp60_1 = B_1_0*w39;                                      const double tmp60_1 = B_1_0*w39;
2708                                      const register double tmp32_1 = tmp1_0*w36;                                      const double tmp32_1 = tmp1_0*w36;
2709                                      const register double tmp10_1 = B_0_1*w37;                                      const double tmp10_1 = B_0_1*w37;
2710                                      const register double tmp14_1 = B_0_0*w35;                                      const double tmp14_1 = B_0_0*w35;
2711                                      const register double tmp29_1 = B_0_1*w35;                                      const double tmp29_1 = B_0_1*w35;
2712                                      const register double tmp26_1 = B_1_2*w36;                                      const double tmp26_1 = B_1_2*w36;
2713                                      const register double tmp30_1 = B_1_3*w43;                                      const double tmp30_1 = B_1_3*w43;
2714                                      const register double tmp70_1 = B_0_2*w35;                                      const double tmp70_1 = B_0_2*w35;
2715                                      const register double tmp34_1 = B_1_3*w39;                                      const double tmp34_1 = B_1_3*w39;
2716                                      const register double tmp51_1 = B_1_0*w43;                                      const double tmp51_1 = B_1_0*w43;
2717                                      const register double tmp31_1 = B_0_2*w34;                                      const double tmp31_1 = B_0_2*w34;
2718                                      const register double tmp45_1 = tmp3_0*w28;                                      const double tmp45_1 = tmp3_0*w28;
2719                                      const register double tmp11_1 = tmp1_0*w39;                                      const double tmp11_1 = tmp1_0*w39;
2720                                      const register double tmp52_1 = B_1_1*w29;                                      const double tmp52_1 = B_1_1*w29;
2721                                      const register double tmp44_1 = B_1_3*w32;                                      const double tmp44_1 = B_1_3*w32;
2722                                      const register double tmp25_1 = B_1_1*w39;                                      const double tmp25_1 = B_1_1*w39;
2723                                      const register double tmp47_1 = tmp2_0*w33;                                      const double tmp47_1 = tmp2_0*w33;
2724                                      const register double tmp72_1 = B_1_3*w29;                                      const double tmp72_1 = B_1_3*w29;
2725                                      const register double tmp40_1 = tmp2_0*w28;                                      const double tmp40_1 = tmp2_0*w28;
2726                                      const register double tmp46_1 = B_1_2*w40;                                      const double tmp46_1 = B_1_2*w40;
2727                                      const register double tmp36_1 = B_1_2*w42;                                      const double tmp36_1 = B_1_2*w42;
2728                                      const register double tmp24_1 = B_0_0*w37;                                      const double tmp24_1 = B_0_0*w37;
2729                                      const register double tmp77_1 = B_0_3*w35;                                      const double tmp77_1 = B_0_3*w35;
2730                                      const register double tmp68_1 = B_0_3*w37;                                      const double tmp68_1 = B_0_3*w37;
2731                                      const register double tmp78_1 = B_0_0*w34;                                      const double tmp78_1 = B_0_0*w34;
2732                                      const register double tmp12_1 = tmp0_0*w36;                                      const double tmp12_1 = tmp0_0*w36;
2733                                      const register double tmp75_1 = B_1_0*w32;                                      const double tmp75_1 = B_1_0*w32;
2734                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2735                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2736                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
# Line 2768  void Rectangle::assemblePDESystem(Paso_S Line 2752  void Rectangle::assemblePDESystem(Paso_S
2752                          } else { // constant data                          } else { // constant data
2753                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2754                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2755                                      const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];                                      const double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
2756                                      const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];                                      const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
2757                                      const register double tmp6_1 = B_1*w50;                                      const double tmp6_1 = B_1*w50;
2758                                      const register double tmp1_1 = B_1*w45;                                      const double tmp1_1 = B_1*w45;
2759                                      const register double tmp5_1 = B_1*w49;                                      const double tmp5_1 = B_1*w49;
2760                                      const register double tmp4_1 = B_1*w48;                                      const double tmp4_1 = B_1*w48;
2761                                      const register double tmp0_1 = B_0*w44;                                      const double tmp0_1 = B_0*w44;
2762                                      const register double tmp2_1 = B_0*w46;                                      const double tmp2_1 = B_0*w46;
2763                                      const register double tmp7_1 = B_0*w51;                                      const double tmp7_1 = B_0*w51;
2764                                      const register double tmp3_1 = B_0*w47;                                      const double tmp3_1 = B_0*w47;
2765                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
2766                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
2767                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
# Line 2807  void Rectangle::assemblePDESystem(Paso_S Line 2791  void Rectangle::assemblePDESystem(Paso_S
2791                          if (C.actsExpanded()) {                          if (C.actsExpanded()) {
2792                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2793                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2794                                      const register double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];                                      const double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
2795                                      const register double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];                                      const double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
2796                                      const register double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];                                      const double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
2797                                      const register double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];                                      const double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
2798                                      const register double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];                                      const double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];
2799                                      const register double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];                                      const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
2800                                      const register double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];                                      const double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];
2801                                      const register double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];                                      const double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];
2802                                      const register double tmp2_0 = C_0_1 + C_0_3;                                      const double tmp2_0 = C_0_1 + C_0_3;
2803                                      const register double tmp1_0 = C_1_2 + C_1_3;                                      const double tmp1_0 = C_1_2 + C_1_3;
2804                                      const register double tmp3_0 = C_0_0 + C_0_2;                                      const double tmp3_0 = C_0_0 + C_0_2;
2805                                      const register double tmp0_0 = C_1_0 + C_1_1;                                      const double tmp0_0 = C_1_0 + C_1_1;
2806                                      const register double tmp64_1 = C_0_2*w30;                                      const double tmp64_1 = C_0_2*w30;
2807                                      const register double tmp14_1 = C_0_2*w28;                                      const double tmp14_1 = C_0_2*w28;
2808                                      const register double tmp19_1 = C_0_3*w31;                                      const double tmp19_1 = C_0_3*w31;
2809                                      const register double tmp22_1 = C_1_0*w40;                                      const double tmp22_1 = C_1_0*w40;
2810                                      const register double tmp37_1 = tmp3_0*w35;                                      const double tmp37_1 = tmp3_0*w35;
2811                                      const register double tmp29_1 = C_0_1*w35;                                      const double tmp29_1 = C_0_1*w35;
2812                                      const register double tmp73_1 = C_0_2*w37;                                      const double tmp73_1 = C_0_2*w37;
2813                                      const register double tmp74_1 = C_1_2*w41;                                      const double tmp74_1 = C_1_2*w41;
2814                                      const register double tmp52_1 = C_1_3*w39;                                      const double tmp52_1 = C_1_3*w39;
2815                                      const register double tmp25_1 = C_1_1*w39;                                      const double tmp25_1 = C_1_1*w39;
2816                                      const register double tmp62_1 = C_0_1*w31;                                      const double tmp62_1 = C_0_1*w31;
2817                                      const register double tmp79_1 = C_1_1*w40;                                      const double tmp79_1 = C_1_1*w40;
2818                                      const register double tmp43_1 = C_1_1*w36;                                      const double tmp43_1 = C_1_1*w36;
2819                                      const register double tmp27_1 = C_0_3*w38;                                      const double tmp27_1 = C_0_3*w38;
2820                                      const register double tmp28_1 = C_1_0*w42;                                      const double tmp28_1 = C_1_0*w42;
2821                                      const register double tmp63_1 = C_1_1*w42;                                      const double tmp63_1 = C_1_1*w42;
2822                                      const register double tmp59_1 = C_0_3*w34;                                      const double tmp59_1 = C_0_3*w34;
2823                                      const register double tmp72_1 = C_1_3*w29;                                      const double tmp72_1 = C_1_3*w29;
2824                                      const register double tmp40_1 = tmp2_0*w35;                                      const double tmp40_1 = tmp2_0*w35;
2825                                      const register double tmp13_1 = C_0_3*w30;                                      const double tmp13_1 = C_0_3*w30;
2826                                      const register double tmp51_1 = C_1_2*w40;                                      const double tmp51_1 = C_1_2*w40;
2827                                      const register double tmp54_1 = C_1_2*w42;                                      const double tmp54_1 = C_1_2*w42;
2828                                      const register double tmp12_1 = C_0_0*w31;                                      const double tmp12_1 = C_0_0*w31;
2829                                      const register double tmp2_1 = tmp1_0*w32;                                      const double tmp2_1 = tmp1_0*w32;
2830                                      const register double tmp68_1 = C_0_2*w31;                                      const double tmp68_1 = C_0_2*w31;
2831                                      const register double tmp75_1 = C_1_0*w32;                                      const double tmp75_1 = C_1_0*w32;
2832                                      const register double tmp49_1 = C_1_1*w41;                                      const double tmp49_1 = C_1_1*w41;
2833                                      const register double tmp4_1 = C_0_2*w35;                                      const double tmp4_1 = C_0_2*w35;
2834                                      const register double tmp66_1 = C_0_3*w28;                                      const double tmp66_1 = C_0_3*w28;
2835                                      const register double tmp56_1 = C_0_1*w37;                                      const double tmp56_1 = C_0_1*w37;
2836                                      const register double tmp5_1 = C_0_1*w34;                                      const double tmp5_1 = C_0_1*w34;
2837                                      const register double tmp38_1 = tmp2_0*w34;                                      const double tmp38_1 = tmp2_0*w34;
2838                                      const register double tmp76_1 = C_0_1*w38;                                      const double tmp76_1 = C_0_1*w38;
2839                                      const register double tmp21_1 = C_0_1*w28;                                      const double tmp21_1 = C_0_1*w28;
2840                                      const register double tmp69_1 = C_0_1*w30;                                      const double tmp69_1 = C_0_1*w30;
2841                                      const register double tmp53_1 = C_1_0*w36;                                      const double tmp53_1 = C_1_0*w36;
2842                                      const register double tmp42_1 = C_1_2*w39;                                      const double tmp42_1 = C_1_2*w39;
2843                                      const register double tmp32_1 = tmp1_0*w29;                                      const double tmp32_1 = tmp1_0*w29;
2844                                      const register double tmp45_1 = C_1_0*w43;                                      const double tmp45_1 = C_1_0*w43;
2845                                      const register double tmp33_1 = tmp0_0*w32;                                      const double tmp33_1 = tmp0_0*w32;
2846                                      const register double tmp35_1 = C_1_0*w41;                                      const double tmp35_1 = C_1_0*w41;
2847                                      const register double tmp26_1 = C_1_2*w36;                                      const double tmp26_1 = C_1_2*w36;
2848                                      const register double tmp67_1 = C_0_0*w33;                                      const double tmp67_1 = C_0_0*w33;
2849                                      const register double tmp31_1 = C_0_2*w34;                                      const double tmp31_1 = C_0_2*w34;
2850                                      const register double tmp20_1 = C_0_0*w30;                                      const double tmp20_1 = C_0_0*w30;
2851                                      const register double tmp70_1 = C_0_0*w28;                                      const double tmp70_1 = C_0_0*w28;
2852                                      const register double tmp8_1 = tmp0_0*w39;                                      const double tmp8_1 = tmp0_0*w39;
2853                                      const register double tmp30_1 = C_1_3*w43;                                      const double tmp30_1 = C_1_3*w43;
2854                                      const register double tmp0_1 = tmp0_0*w29;                                      const double tmp0_1 = tmp0_0*w29;
2855                                      const register double tmp17_1 = C_1_3*w41;                                      const double tmp17_1 = C_1_3*w41;
2856                                      const register double tmp58_1 = C_0_0*w35;                                      const double tmp58_1 = C_0_0*w35;
2857                                      const register double tmp9_1 = tmp3_0*w33;                                      const double tmp9_1 = tmp3_0*w33;
2858                                      const register double tmp61_1 = C_1_3*w36;                                      const double tmp61_1 = C_1_3*w36;
2859                                      const register double tmp41_1 = tmp3_0*w34;                                      const double tmp41_1 = tmp3_0*w34;
2860                                      const register double tmp50_1 = C_1_3*w32;                                      const double tmp50_1 = C_1_3*w32;
2861                                      const register double tmp18_1 = C_1_1*w32;                                      const double tmp18_1 = C_1_1*w32;
2862                                      const register double tmp6_1 = tmp1_0*w36;                                      const double tmp6_1 = tmp1_0*w36;
2863                                      const register double tmp3_1 = C_0_0*w38;                                      const double tmp3_1 = C_0_0*w38;
2864                                      const register double tmp34_1 = C_1_1*w29;                                      const double tmp34_1 = C_1_1*w29;
2865                                      const register double tmp77_1 = C_0_3*w35;                                      const double tmp77_1 = C_0_3*w35;
2866                                      const register double tmp65_1 = C_1_2*w43;                                      const double tmp65_1 = C_1_2*w43;
2867                                      const register double tmp71_1 = C_0_3*w33;                                      const double tmp71_1 = C_0_3*w33;
2868                                      const register double tmp55_1 = C_1_1*w43;                                      const double tmp55_1 = C_1_1*w43;
2869                                      const register double tmp46_1 = tmp3_0*w28;                                      const double tmp46_1 = tmp3_0*w28;
2870                                      const register double tmp24_1 = C_0_0*w37;                                      const double tmp24_1 = C_0_0*w37;
2871                                      const register double tmp10_1 = tmp1_0*w39;                                      const double tmp10_1 = tmp1_0*w39;
2872                                      const register double tmp48_1 = C_1_0*w29;                                      const double tmp48_1 = C_1_0*w29;
2873                                      const register double tmp15_1 = C_0_1*w33;                                      const double tmp15_1 = C_0_1*w33;
2874                                      const register double tmp36_1 = C_1_2*w32;                                      const double tmp36_1 = C_1_2*w32;
2875                                      const register double tmp60_1 = C_1_0*w39;                                      const double tmp60_1 = C_1_0*w39;
2876                                      const register double tmp47_1 = tmp2_0*w33;                                      const double tmp47_1 = tmp2_0*w33;
2877                                      const register double tmp16_1 = C_1_2*w29;                                      const double tmp16_1 = C_1_2*w29;
2878                                      const register double tmp1_1 = C_0_3*w37;                                      const double tmp1_1 = C_0_3*w37;
2879                                      const register double tmp7_1 = tmp2_0*w28;                                      const double tmp7_1 = tmp2_0*w28;
2880                                      const register double tmp39_1 = C_1_3*w40;                                      const double tmp39_1 = C_1_3*w40;
2881                                      const register double tmp44_1 = C_1_3*w42;                                      const double tmp44_1 = C_1_3*w42;
2882                                      const register double tmp57_1 = C_0_2*w38;                                      const double tmp57_1 = C_0_2*w38;
2883                                      const register double tmp78_1 = C_0_0*w34;                                      const double tmp78_1 = C_0_0*w34;
2884                                      const register double tmp11_1 = tmp0_0*w36;                                      const double tmp11_1 = tmp0_0*w36;
2885                                      const register double tmp23_1 = C_0_2*w33;                                      const double tmp23_1 = C_0_2*w33;
2886                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2887                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2888                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
# Line 2920  void Rectangle::assemblePDESystem(Paso_S Line 2904  void Rectangle::assemblePDESystem(Paso_S
2904                          } else { // constant data                          } else { // constant data
2905                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2906                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2907                                      const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
2908                                      const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
2909                                      const register double tmp1_1 = C_1*w45;                                      const double tmp1_1 = C_1*w45;
2910                                      const register double tmp3_1 = C_0*w51;                                      const double tmp3_1 = C_0*w51;
2911                                      const register double tmp4_1 = C_0*w44;                                      const double tmp4_1 = C_0*w44;
2912                                      const register double tmp7_1 = C_0*w46;                                      const double tmp7_1 = C_0*w46;
2913                                      const register double tmp5_1 = C_1*w49;                                      const double tmp5_1 = C_1*w49;
2914                                      const register double tmp2_1 = C_1*w48;                                      const double tmp2_1 = C_1*w48;
2915                                      const register double tmp0_1 = C_0*w47;                                      const double tmp0_1 = C_0*w47;
2916                                      const register double tmp6_1 = C_1*w50;                                      const double tmp6_1 = C_1*w50;
2917                                      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;
2918                                      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;
2919                                      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 2943  void Rectangle::assemblePDESystem(Paso_S
2943                          if (D.actsExpanded()) {                          if (D.actsExpanded()) {
2944                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2945                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2946                                      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)];
2947                                      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)];
2948                                      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)];
2949                                      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)];
2950                                      const register double tmp4_0 = D_1 + D_3;                                      const double tmp4_0 = D_1 + D_3;
2951                                      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;
2952                                      const register double tmp5_0 = D_0 + D_2;                                      const double tmp5_0 = D_0 + D_2;
2953                                      const register double tmp0_0 = D_0 + D_1;                                      const double tmp0_0 = D_0 + D_1;
2954                                      const register double tmp6_0 = D_0 + D_3;                                      const double tmp6_0 = D_0 + D_3;
2955                                      const register double tmp1_0 = D_2 + D_3;                                      const double tmp1_0 = D_2 + D_3;
2956                                      const register double tmp3_0 = D_1 + D_2;                                      const double tmp3_0 = D_1 + D_2;
2957                                      const register double tmp16_1 = D_1*w56;                                      const double tmp16_1 = D_1*w56;
2958                                      const register double tmp14_1 = tmp6_0*w54;                                      const double tmp14_1 = tmp6_0*w54;
2959                                      const register double tmp8_1 = D_3*w55;                                      const double tmp8_1 = D_3*w55;
2960                                      const register double tmp2_1 = tmp2_0*w54;                                      const double tmp2_1 = tmp2_0*w54;
2961                                      const register double tmp12_1 = tmp5_0*w52;                                      const double tmp12_1 = tmp5_0*w52;
2962                                      const register double tmp4_1 = tmp0_0*w53;                                      const double tmp4_1 = tmp0_0*w53;
2963                                      const register double tmp3_1 = tmp1_0*w52;                                      const double tmp3_1 = tmp1_0*w52;
2964                                      const register double tmp13_1 = tmp4_0*w53;                                      const double tmp13_1 = tmp4_0*w53;
2965                                      const register double tmp10_1 = tmp4_0*w52;                                      const double tmp10_1 = tmp4_0*w52;
2966                                      const register double tmp1_1 = tmp1_0*w53;                                      const double tmp1_1 = tmp1_0*w53;
2967                                      const register double tmp7_1 = D_3*w56;                                      const double tmp7_1 = D_3*w56;
2968                                      const register double tmp0_1 = tmp0_0*w52;                                      const double tmp0_1 = tmp0_0*w52;
2969                                      const register double tmp11_1 = tmp5_0*w53;                                      const double tmp11_1 = tmp5_0*w53;
2970                                      const register double tmp9_1 = D_0*w56;                                      const double tmp9_1 = D_0*w56;
2971                                      const register double tmp5_1 = tmp3_0*w54;                                      const double tmp5_1 = tmp3_0*w54;
2972                                      const register double tmp18_1 = D_2*w56;                                      const double tmp18_1 = D_2*w56;
2973                                      const register double tmp17_1 = D_1*w55;                                      const double tmp17_1 = D_1*w55;
2974                                      const register double tmp6_1 = D_0*w55;                                      const double tmp6_1 = D_0*w55;
2975                                      const register double tmp15_1 = D_2*w55;                                      const double tmp15_1 = D_2*w55;
2976                                      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;
2977                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1;
2978                                      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 2994  void Rectangle::assemblePDESystem(Paso_S
2994                          } else { // constant data                          } else { // constant data
2995                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2996                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2997                                      const register double D_0 = D_p[INDEX2(k, m, numEq)];                                      const double D_0 = D_p[INDEX2(k, m, numEq)];
2998                                      const register double tmp2_1 = D_0*w59;                                      const double tmp2_1 = D_0*w59;
2999                                      const register double tmp1_1 = D_0*w58;                                      const double tmp1_1 = D_0*w58;
3000                                      const register double tmp0_1 = D_0*w57;                                      const double tmp0_1 = D_0*w57;
3001                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
3002                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1;
3003                                      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 3026  void Rectangle::assemblePDESystem(Paso_S
3026                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3027                          if (X.actsExpanded()) {                          if (X.actsExpanded()) {
3028                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3029                                  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)];
3030                                  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)];
3031                                  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)];
3032                                  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)];
3033                                  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)];
3034                                  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)];
3035                                  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)];
3036                                  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)];
3037                                  const register double tmp1_0 = X_1_1 + X_1_3;                                  const double tmp1_0 = X_1_1 + X_1_3;
3038                                  const register double tmp3_0 = X_0_0 + X_0_1;                                  const double tmp3_0 = X_0_0 + X_0_1;
3039                                  const register double tmp2_0 = X_1_0 + X_1_2;                                  const double tmp2_0 = X_1_0 + X_1_2;
3040                                  const register double tmp0_0 = X_0_2 + X_0_3;                                  const double tmp0_0 = X_0_2 + X_0_3;
3041                                  const register double tmp8_1 = tmp2_0*w66;                                  const double tmp8_1 = tmp2_0*w66;
3042                                  const register double tmp5_1 = tmp3_0*w64;                                  const double tmp5_1 = tmp3_0*w64;
3043                                  const register double tmp14_1 = tmp0_0*w64;                                  const double tmp14_1 = tmp0_0*w64;
3044                                  const register double tmp3_1 = tmp3_0*w60;                                  const double tmp3_1 = tmp3_0*w60;
3045                                  const register double tmp9_1 = tmp3_0*w63;                                  const double tmp9_1 = tmp3_0*w63;
3046                                  const register double tmp13_1 = tmp3_0*w65;                                  const double tmp13_1 = tmp3_0*w65;
3047                                  const register double tmp12_1 = tmp1_0*w66;                                  const double tmp12_1 = tmp1_0*w66;
3048                                  const register double tmp10_1 = tmp0_0*w60;                                  const double tmp10_1 = tmp0_0*w60;
3049                                  const register double tmp2_1 = tmp2_0*w61;                                  const double tmp2_1 = tmp2_0*w61;
3050                                  const register double tmp6_1 = tmp2_0*w62;                                  const double tmp6_1 = tmp2_0*w62;
3051                                  const register double tmp4_1 = tmp0_0*w65;                                  const double tmp4_1 = tmp0_0*w65;
3052                                  const register double tmp11_1 = tmp1_0*w67;                                  const double tmp11_1 = tmp1_0*w67;
3053                                  const register double tmp1_1 = tmp1_0*w62;                                  const double tmp1_1 = tmp1_0*w62;
3054                                  const register double tmp7_1 = tmp1_0*w61;                                  const double tmp7_1 = tmp1_0*w61;
3055                                  const register double tmp0_1 = tmp0_0*w63;                                  const double tmp0_1 = tmp0_0*w63;
3056                                  const register double tmp15_1 = tmp2_0*w67;                                  const double tmp15_1 = tmp2_0*w67;
3057                                  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;
3058                                  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;
3059                                  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 3061  void Rectangle::assemblePDESystem(Paso_S
3061                              }                              }
3062                          } else { // constant data                          } else { // constant data
3063                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3064                                  const register double X_0 = X_p[INDEX2(k, 0, numEq)];                                  const double X_0 = X_p[INDEX2(k, 0, numEq)];
3065                                  const register double X_1 = X_p[INDEX2(k, 1, numEq)];                                  const double X_1 = X_p[INDEX2(k, 1, numEq)];
3066                                  const register double tmp0_1 = X_1*w69;                                  const double tmp0_1 = X_1*w69;
3067                                  const register double tmp1_1 = X_0*w68;                                  const double tmp1_1 = X_0*w68;
3068                                  const register double tmp2_1 = X_0*w70;                                  const double tmp2_1 = X_0*w70;
3069                                  const register double tmp3_1 = X_1*w71;                                  const double tmp3_1 = X_1*w71;
3070                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3071                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;
3072                                  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 3082  void Rectangle::assemblePDESystem(Paso_S
3082                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3083                          if (Y.actsExpanded()) {                          if (Y.actsExpanded()) {
3084                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3085                                  const register double Y_0 = Y_p[INDEX2(k, 0, numEq)];                                  const double Y_0 = Y_p[INDEX2(k, 0, numEq)];
3086                                  const register double Y_1 = Y_p[INDEX2(k, 1, numEq)];                                  const double Y_1 = Y_p[INDEX2(k, 1, numEq)];
3087                                  const register double Y_2 = Y_p[INDEX2(k, 2, numEq)];                                  const double Y_2 = Y_p[INDEX2(k, 2, numEq)];
3088                                  const register double Y_3 = Y_p[INDEX2(k, 3, numEq)];                                  const double Y_3 = Y_p[INDEX2(k, 3, numEq)];
3089                                  const register double tmp0_0 = Y_1 + Y_2;                                  const double tmp0_0 = Y_1 + Y_2;
3090                                  const register double tmp1_0 = Y_0 + Y_3;                                  const double tmp1_0 = Y_0 + Y_3;
3091                                  const register double tmp0_1 = Y_0*w72;                                  const double tmp0_1 = Y_0*w72;
3092                                  const register double tmp1_1 = tmp0_0*w73;                                  const double tmp1_1 = tmp0_0*w73;
3093                                  const register double tmp2_1 = Y_3*w74;                                  const double tmp2_1 = Y_3*w74;
3094                                  const register double tmp3_1 = Y_1*w72;                                  const double tmp3_1 = Y_1*w72;
3095                                  const register double tmp4_1 = tmp1_0*w73;                                  const double tmp4_1 = tmp1_0*w73;
3096                                  const register double tmp5_1 = Y_2*w74;                                  const double tmp5_1 = Y_2*w74;
3097                                  const register double tmp6_1 = Y_2*w72;                                  const double tmp6_1 = Y_2*w72;
3098                                  const register double tmp7_1 = Y_1*w74;                                  const double tmp7_1 = Y_1*w74;
3099                                  const register double tmp8_1 = Y_3*w72;                                  const double tmp8_1 = Y_3*w72;
3100                                  const register double tmp9_1 = Y_0*w74;                                  const double tmp9_1 = Y_0*w74;
3101                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
3102                                  EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
3103                                  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 3105  void Rectangle::assemblePDESystem(Paso_S
3105                              }                              }
3106                          } else { // constant data                          } else { // constant data
3107                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3108                                  const register double tmp0_1 = Y_p[k]*w75;                                  const double tmp0_1 = Y_p[k]*w75;
3109                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3110                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3111                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
# Line 3161  void Rectangle::assemblePDESystem(Paso_S Line 3145  void Rectangle::assemblePDESystem(Paso_S
3145  void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,
3146          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
3147          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
3148          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  
3149  {  {
3150      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
3151      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 3211  void Rectangle::assemblePDESystemReduced Line 3194  void Rectangle::assemblePDESystemReduced
3194                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
3195                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3196                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3197                                  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)];
3198                                  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)];
3199                                  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)];
3200                                  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)];
3201                                  const register double tmp0_0 = A_01 + A_10;                                  const double tmp0_0 = A_01 + A_10;
3202                                  const register double tmp0_1 = A_11*w3;                                  const double tmp0_1 = A_11*w3;
3203                                  const register double tmp1_1 = A_00*w0;                                  const double tmp1_1 = A_00*w0;
3204                                  const register double tmp2_1 = A_01*w1;                                  const double tmp2_1 = A_01*w1;
3205                                  const register double tmp3_1 = A_10*w2;                                  const double tmp3_1 = A_10*w2;
3206                                  const register double tmp4_1 = tmp0_0*w1;                                  const double tmp4_1 = tmp0_0*w1;
3207                                  const register double tmp5_1 = A_11*w4;                                  const double tmp5_1 = A_11*w4;
3208                                  const register double tmp6_1 = A_00*w5;                                  const double tmp6_1 = A_00*w5;
3209                                  const register double tmp7_1 = tmp0_0*w2;                                  const double tmp7_1 = tmp0_0*w2;
3210                                  const register double tmp8_1 = A_10*w1;                                  const double tmp8_1 = A_10*w1;
3211                                  const register double tmp9_1 = A_01*w2;                                  const double tmp9_1 = A_01*w2;
3212                                  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;
3213                                  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;
3214                                  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 3236  void Rectangle::assemblePDESystemReduced
3236                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
3237                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3238                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3239                                  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)];
3240                                  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)];
3241                                  const register double tmp0_1 = B_0*w6;                                  const double tmp0_1 = B_0*w6;
3242                                  const register double tmp1_1 = B_1*w7;                                  const double tmp1_1 = B_1*w7;
3243                                  const register double tmp2_1 = B_0*w8;                                  const double tmp2_1 = B_0*w8;
3244                                  const register double tmp3_1 = B_1*w9;                                  const double tmp3_1 = B_1*w9;
3245                                  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;
3246                                  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;
3247                                  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 3269  void Rectangle::assemblePDESystemReduced
3269                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
3270                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3271                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3272                                  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)];
3273                                  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)];
3274                                  const register double tmp0_1 = C_0*w8;                                  const double tmp0_1 = C_0*w8;
3275                                  const register double tmp1_1 = C_1*w7;                                  const double tmp1_1 = C_1*w7;
3276                                  const register double tmp2_1 = C_1*w9;                                  const double tmp2_1 = C_1*w9;
3277                                  const register double tmp3_1 = C_0*w6;                                  const double tmp3_1 = C_0*w6;
3278                                  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;
3279                                  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;
3280                                  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 3302  void Rectangle::assemblePDESystemReduced
3302                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
3303                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3304                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3305                                  const register double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;                                  const double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;
3306                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
3307                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
3308                                  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 3329  void Rectangle::assemblePDESystemReduced
3329                          add_EM_F=true;                          add_EM_F=true;
3330                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3331                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3332                              const register double X_0 = X_p[INDEX2(k, 0, numEq)];                              const double X_0 = X_p[INDEX2(k, 0, numEq)];
3333                              const register double X_1 = X_p[INDEX2(k, 1, numEq)];                              const double X_1 = X_p[INDEX2(k, 1, numEq)];
3334                              const register double tmp0_1 = X_0*w11;                              const double tmp0_1 = X_0*w11;
3335                              const register double tmp1_1 = X_1*w12;                              const double tmp1_1 = X_1*w12;
3336                              const register double tmp2_1 = X_0*w13;                              const double tmp2_1 = X_0*w13;
3337                              const register double tmp3_1 = X_1*w14;                              const double tmp3_1 = X_1*w14;
3338                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3339                              EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;
3340                              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 3348  void Rectangle::assemblePDESystemReduced
3348                          add_EM_F=true;                          add_EM_F=true;
3349                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3350                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3351                              const register double tmp0_1 = Y_p[k]*w15;                              const double tmp0_1 = Y_p[k]*w15;
3352                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3353                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3354                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
# Line 3400  void Rectangle::assemblePDESystemReduced Line 3383  void Rectangle::assemblePDESystemReduced
3383      } // end of parallel region      } // end of parallel region
3384  }  }
3385    
3386    //protected
3387    void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
3388          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3389    {
3390    }
3391    
3392    //protected
3393    void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
3394          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3395    {
3396    }
3397    
3398    //protected
3399    void Rectangle::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
3400          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3401    {
3402    }
3403    
3404    //protected
3405    void Rectangle::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
3406          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3407    {
3408    }
3409    
3410    
3411  } // end of namespace ripley  } // end of namespace ripley
3412    

Legend:
Removed from v.3768  
changed lines
  Added in v.3769

  ViewVC Help
Powered by ViewVC 1.1.26