/[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 3733 by caltinay, Fri Dec 9 04:02:56 2011 UTC revision 3756 by caltinay, Fri Jan 6 02:35:19 2012 UTC
# Line 13  Line 13 
13    
14  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
15  extern "C" {  extern "C" {
16  #include "paso/SystemMatrixPattern.h"  #include "paso/SystemMatrix.h"
17  }  }
18    
19  #if USE_SILO  #if USE_SILO
# Line 43  Rectangle::Rectangle(int n0, int n1, dou Line 43  Rectangle::Rectangle(int n0, int n1, dou
43      if (m_NX*m_NY != m_mpiInfo->size)      if (m_NX*m_NY != m_mpiInfo->size)
44          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
45    
46      if (n0%m_NX > 0 || n1%m_NY > 0)      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)
47          throw RipleyException("Number of elements must be separable into number of ranks in each dimension");          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");
48    
49      // local number of elements      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
50      m_NE0 = n0/m_NX;          throw RipleyException("Too few elements for the number of ranks");
51      m_NE1 = n1/m_NY;  
52      // local number of nodes (not necessarily owned)      // local number of elements (including overlap)
53        m_NE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
54        if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
55            m_NE0++;
56        m_NE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
57        if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)
58            m_NE1++;
59    
60        // local number of nodes
61      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
62      m_N1 = m_NE1+1;      m_N1 = m_NE1+1;
63    
64      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
65      m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);
66      m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);      if (m_offset0 > 0)
67            m_offset0--;
68        m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);
69        if (m_offset1 > 0)
70            m_offset1--;
71    
72      populateSampleIds();      populateSampleIds();
73        createPattern();
74  }  }
75    
76  Rectangle::~Rectangle()  Rectangle::~Rectangle()
# Line 69  string Rectangle::getDescription() const Line 84  string Rectangle::getDescription() const
84    
85  bool Rectangle::operator==(const AbstractDomain& other) const  bool Rectangle::operator==(const AbstractDomain& other) const
86  {  {
87      if (dynamic_cast<const Rectangle*>(&other))      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
88          return this==&other;      if (o) {
89            return (RipleyDomain::operator==(other) &&
90                    m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1
91                    && m_l0==o->m_l0 && m_l1==o->m_l1
92                    && m_NX==o->m_NX && m_NY==o->m_NY);
93        }
94    
95      return false;      return false;
96  }  }
# Line 221  const int* Rectangle::borrowSampleRefere Line 241  const int* Rectangle::borrowSampleRefere
241  {  {
242      switch (fsType) {      switch (fsType) {
243          case Nodes:          case Nodes:
244          case DegreesOfFreedom: //FIXME          case ReducedNodes: //FIXME: reduced
         case ReducedDegreesOfFreedom: //FIXME  
245              return &m_nodeId[0];              return &m_nodeId[0];
246            case DegreesOfFreedom:
247            case ReducedDegreesOfFreedom: //FIXME: reduced
248                return &m_dofId[0];
249          case Elements:          case Elements:
250          case ReducedElements:          case ReducedElements:
251              return &m_elementId[0];              return &m_elementId[0];
# Line 243  const int* Rectangle::borrowSampleRefere Line 265  const int* Rectangle::borrowSampleRefere
265  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsCode, index_t id) const
266  {  {
267  #ifdef ESYS_MPI  #ifdef ESYS_MPI
268      if (fsCode != ReducedNodes) {      if (fsCode == Nodes) {
269          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];          return (m_dofMap[id] < getNumDOF());
         const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;  
         return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);  
270      } else {      } else {
271          stringstream msg;          stringstream msg;
272          msg << "ownSample() not implemented for "          msg << "ownSample() not implemented for "
# Line 724  void Rectangle::setToNormal(escript::Dat Line 744  void Rectangle::setToNormal(escript::Dat
744      }      }
745  }  }
746    
 void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,  
         escript::Data& rhs, const escript::Data& A, const escript::Data& B,  
         const escript::Data& C, const escript::Data& D,  
         const escript::Data& X, const escript::Data& Y,  
         const escript::Data& d, const escript::Data& y,  
         const escript::Data& d_contact, const escript::Data& y_contact,  
         const escript::Data& d_dirac, const escript::Data& y_dirac) const  
 {  
     throw RipleyException("addPDEToSystem() not implemented");  
 }  
   
747  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
748                                                  bool reducedColOrder) const                                                  bool reducedColOrder) const
749  {  {
750      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
751          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
752    
753      // connector      return m_pattern;
     RankVector neighbour;  
     IndexVector offsetInShared(1,0);  
     IndexVector sendShared, recvShared;  
     const IndexVector faces=getNumFacesPerBoundary();  
     const index_t left = (m_offset0==0 ? 0 : 1);  
     const index_t bottom = (m_offset1==0 ? 0 : 1);  
     // corner node from bottom-left  
     if (faces[0] == 0 && faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX-1);  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0+left]);  
         recvShared.push_back(m_nodeId[0]);  
     }  
     // bottom edge  
     if (faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX);  
         offsetInShared.push_back(offsetInShared.back()+m_N0-left);  
         for (dim_t i=left; i<m_N0; i++) {  
             // easy case, we know the neighbour id's  
             sendShared.push_back(m_nodeId[m_N0+i]);  
             recvShared.push_back(m_nodeId[i]);  
         }  
     }  
     // corner node from bottom-right  
     if (faces[1] == 0 && faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX+1);  
         const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour.back()/m_NX == 0 ? m_N1 : m_N1-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[(bottom+1)*m_N0-1]);  
         recvShared.push_back(first+N0*(N1-1));  
     }  
     // left edge  
     if (faces[0] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-1);  
         offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);  
         for (dim_t i=bottom; i<m_N1; i++) {  
             // easy case, we know the neighbour id's  
             sendShared.push_back(m_nodeId[i*m_N0+1]);  
             recvShared.push_back(m_nodeId[i*m_N0]);  
         }  
     }  
     // right edge  
     if (faces[1] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+1);  
         const index_t rightN0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);  
         for (dim_t i=bottom; i<m_N1; i++) {  
             sendShared.push_back(m_nodeId[(i+1)*m_N0-1]);  
             recvShared.push_back(first+rightN0*(i-bottom));  
         }  
     }  
     // corner node from top-left  
     if (faces[0] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX-1);  
         const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+left]);  
         recvShared.push_back(first+N0-1);  
     }  
     // top edge  
     if (faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+m_N0-left);  
         for (dim_t i=left; i<m_N0; i++) {  
             sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+i]);  
             recvShared.push_back(first+i-left);  
         }  
     }  
     // corner node from top-right  
     if (faces[1] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX+1);  
         const index_t first=m_nodeDistribution[neighbour.back()];  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(m_nodeId[m_N0*m_N1-1]);  
         recvShared.push_back(first);  
     }  
     const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];  
     /*  
     cout << "--- rcv_shcomp ---" << endl;  
     cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;  
     for (size_t i=0; i<neighbour.size(); i++) {  
         cout << "neighbor[" << i << "]=" << neighbour[i]  
             << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;  
     }  
     for (size_t i=0; i<recvShared.size(); i++) {  
         cout << "shared[" << i << "]=" << recvShared[i] << endl;  
     }  
     cout << "--- snd_shcomp ---" << endl;  
     for (size_t i=0; i<sendShared.size(); i++) {  
         cout << "shared[" << i << "]=" << sendShared[i] << endl;  
     }  
     */  
   
     Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &sendShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &recvShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);  
     Paso_SharedComponents_free(snd_shcomp);  
     Paso_SharedComponents_free(rcv_shcomp);  
   
     // create patterns  
     dim_t M, N;  
     IndexVector ptr(1,0);  
     IndexVector index;  
   
     // main pattern  
     for (index_t i=0; i<numDOF; i++) {  
         // always add the node itself  
         index.push_back(i);  
         int num=insertNeighbours(index, i);  
         ptr.push_back(ptr.back()+num+1);  
     }  
     M=N=ptr.size()-1;  
     // paso will manage the memory  
     index_t* indexC = MEMALLOC(index.size(),index_t);  
     index_t* ptrC = MEMALLOC(ptr.size(), index_t);  
     copy(index.begin(), index.end(), indexC);  
     copy(ptr.begin(), ptr.end(), ptrC);  
     Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
             M, N, ptrC, indexC);  
   
     /*  
     cout << "--- main_pattern ---" << endl;  
     cout << "M=" << M << ", N=" << N << endl;  
     for (size_t i=0; i<ptr.size(); i++) {  
         cout << "ptr[" << i << "]=" << ptr[i] << endl;  
     }  
     for (size_t i=0; i<index.size(); i++) {  
         cout << "index[" << i << "]=" << index[i] << endl;  
     }  
     */  
   
     ptr.clear();  
     index.clear();  
   
     // column & row couple patterns  
     Paso_Pattern *colCouplePattern, *rowCouplePattern;  
     generateCouplePatterns(&colCouplePattern, &rowCouplePattern);  
   
     // allocate paso distribution  
     Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,  
             const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);  
   
     Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(  
             MATRIX_FORMAT_DEFAULT, distribution, distribution,  
             mainPattern, colCouplePattern, rowCouplePattern,  
             connector, connector);  
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colCouplePattern);  
     Paso_Pattern_free(rowCouplePattern);  
     Paso_Distribution_free(distribution);  
     return pattern;  
754  }  }
755    
756  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
# Line 966  pair<double,double> Rectangle::getFirstC Line 815  pair<double,double> Rectangle::getFirstC
815  }  }
816    
817  //protected  //protected
818    dim_t Rectangle::getNumDOF() const
819    {
820        return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;
821    }
822    
823    //protected
824  dim_t Rectangle::getNumFaceElements() const  dim_t Rectangle::getNumFaceElements() const
825  {  {
826      const IndexVector faces = getNumFacesPerBoundary();      const IndexVector faces = getNumFacesPerBoundary();
# Line 998  void Rectangle::assembleCoordinates(escr Line 853  void Rectangle::assembleCoordinates(escr
853      }      }
854  }  }
855    
856  //private  //protected
857  void Rectangle::populateSampleIds()  dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
858  {  {
859      // identifiers are ordered from left to right, bottom to top on each rank,      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
860      // except for the shared nodes which are owned by the rank below / to the      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
861      // left of the current rank      const int x=node%nDOF0;
862        const int y=node/nDOF0;
863      // build node distribution vector first.      dim_t num=0;
864      // m_nodeDistribution[i] is the first node id on rank i, that is      // loop through potential neighbours and add to index if positions are
865      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // within bounds
866      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      for (int i1=-1; i1<2; i1++) {
867      m_nodeDistribution[1]=getNumNodes();          for (int i0=-1; i0<2; i0++) {
868      for (dim_t k=1; k<m_mpiInfo->size-1; k++) {              // skip node itself
869          const index_t x=k%m_NX;              if (i0==0 && i1==0)
870          const index_t y=k/m_NX;                  continue;
871          index_t numNodes=getNumNodes();              // location of neighbour node
872          if (x>0)              const int nx=x+i0;
873              numNodes-=m_N1;              const int ny=y+i1;
874          if (y>0)              if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {
875              numNodes-=m_N0;                  index.push_back(ny*nDOF0+nx);
876          if (x>0 && y>0)                  num++;
877              numNodes++; // subtracted corner twice -> fix that              }
878          m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;          }
879      }      }
     m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();  
880    
881      m_nodeId.resize(getNumNodes());      return num;
882    }
883    
884    //protected
885    void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const
886    {
887        const dim_t numComp = in.getDataPointSize();
888        out.requireWrite();
889    
     // the bottom row and left column are not owned by this rank so the  
     // identifiers need to be computed accordingly  
890      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
891      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
892      if (left>0) {      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
893          const int neighbour=m_mpiInfo->rank-1;      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
         const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
894  #pragma omp parallel for  #pragma omp parallel for
895          for (dim_t i1=bottom; i1<m_N1; i1++) {      for (index_t i=0; i<nDOF1; i++) {
896              m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]          for (index_t j=0; j<nDOF0; j++) {
897                  + (i1-bottom+1)*leftN0-1;              const index_t n=j+left+(i+bottom)*m_N0;
898                const double* src=in.getSampleDataRO(n);
899                copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
900          }          }
901      }      }
902      if (bottom>0) {  }
903          const int neighbour=m_mpiInfo->rank-m_NX;  
904          const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  //protected
905          const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);  void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const
906    {
907        const dim_t numComp = in.getDataPointSize();
908        Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
909        in.requireWrite();
910        Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));
911    
912        const dim_t numDOF = getNumDOF();
913        out.requireWrite();
914        const double* buffer = Paso_Coupler_finishCollect(coupler);
915    
916  #pragma omp parallel for  #pragma omp parallel for
917          for (dim_t i0=left; i0<m_N0; i0++) {      for (index_t i=0; i<getNumNodes(); i++) {
918              m_nodeId[i0]=m_nodeDistribution[neighbour]          const double* src=(m_dofMap[i]<numDOF ?
919                  + (bottomN1-1)*bottomN0 + i0 - left;                  in.getSampleDataRO(m_dofMap[i])
920          }                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
921            copy(src, src+numComp, out.getSampleDataRW(i));
922      }      }
923      if (left>0 && bottom>0) {  }
924          const int neighbour=m_mpiInfo->rank-m_NX-1;  
925          const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  //private
926          const index_t N1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);  void Rectangle::populateSampleIds()
927          m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1-1;  {
928        // identifiers are ordered from left to right, bottom to top globablly.
929    
930        // build node distribution vector first.
931        // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
932        m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
933        const dim_t numDOF=getNumDOF();
934        for (dim_t k=1; k<m_mpiInfo->size; k++) {
935            m_nodeDistribution[k]=k*numDOF;
936      }      }
937        m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
938        m_nodeId.resize(getNumNodes());
939        m_dofId.resize(numDOF);
940        m_elementId.resize(getNumElements());
941        m_faceId.resize(getNumFaceElements());
942    
943      // the rest of the id's are contiguous  #pragma omp parallel
944      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];      {
945  #pragma omp parallel for          // nodes
946      for (dim_t i1=bottom; i1<m_N1; i1++) {  #pragma omp for nowait
947          for (dim_t i0=left; i0<m_N0; i0++) {          for (dim_t i1=0; i1<m_N1; i1++) {
948              m_nodeId[i0+i1*m_N0] = firstId+i0-left+(i1-bottom)*(m_N0-left);              for (dim_t i0=0; i0<m_N0; i0++) {
949                    m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;
950                }
951          }          }
     }  
952    
953      // element id's          // degrees of freedom
954      m_elementId.resize(getNumElements());  #pragma omp for nowait
955  #pragma omp parallel for          for (dim_t k=0; k<numDOF; k++)
956      for (dim_t k=0; k<getNumElements(); k++) {              m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;
         m_elementId[k]=k;  
     }  
957    
958      // face element id's          // elements
959      m_faceId.resize(getNumFaceElements());  #pragma omp for nowait
960  #pragma omp parallel for          for (dim_t i1=0; i1<m_NE1; i1++) {
961      for (dim_t k=0; k<getNumFaceElements(); k++) {              for (dim_t i0=0; i0<m_NE0; i0++) {
962          m_faceId[k]=k;                  m_elementId[i0+i1*m_NE0]=(m_offset1+i1)*m_gNE0+m_offset0+i0;
963      }              }
964            }
965    
966            // face elements
967    #pragma omp for
968            for (dim_t k=0; k<getNumFaceElements(); k++)
969                m_faceId[k]=k;
970        } // end parallel section
971    
972        m_nodeTags.assign(getNumNodes(), 0);
973        updateTagsInUse(Nodes);
974    
975        m_elementTags.assign(getNumElements(), 0);
976        updateTagsInUse(Elements);
977    
978      // generate face offset vector and set face tags      // generate face offset vector and set face tags
979      const IndexVector facesPerEdge = getNumFacesPerBoundary();      const IndexVector facesPerEdge = getNumFacesPerBoundary();
# Line 1101  void Rectangle::populateSampleIds() Line 997  void Rectangle::populateSampleIds()
997  }  }
998    
999  //private  //private
1000  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const  void Rectangle::createPattern()
1001  {  {
1002      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1003      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1004      const int x=node%myN0;      const index_t left = (m_offset0==0 ? 0 : 1);
1005      const int y=node/myN0;      const index_t bottom = (m_offset1==0 ? 0 : 1);
1006      int num=0;  
1007      if (y>0) {      // populate node->DOF mapping with own degrees of freedom.
1008          if (x>0) {      // The rest is assigned in the loop further down
1009              // bottom-left      m_dofMap.assign(getNumNodes(), 0);
1010              index.push_back(node-myN0-1);  #pragma omp parallel for
1011              num++;      for (index_t i=bottom; i<m_N1; i++) {
1012          }          for (index_t j=left; j<m_N0; j++) {
1013          // bottom              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
         index.push_back(node-myN0);  
         num++;  
         if (x<myN0-1) {  
             // bottom-right  
             index.push_back(node-myN0+1);  
             num++;  
         }  
     }  
     if (x<myN0-1) {  
         // right  
         index.push_back(node+1);  
         num++;  
         if (y<myN1-1) {  
             // top-right  
             index.push_back(node+myN0+1);  
             num++;  
1014          }          }
1015      }      }
1016      if (y<myN1-1) {  
1017          // top      // build list of shared components and neighbours by looping through
1018          index.push_back(node+myN0);      // all potential neighbouring ranks and checking if positions are
1019          num++;      // within bounds
1020          if (x>0) {      const dim_t numDOF=nDOF0*nDOF1;
1021              // top-left      vector<IndexVector> colIndices(numDOF); // for the couple blocks
1022              index.push_back(node+myN0-1);      RankVector neighbour;
1023              num++;      IndexVector offsetInShared(1,0);
1024        IndexVector sendShared, recvShared;
1025        int numShared=0;
1026        const int x=m_mpiInfo->rank%m_NX;
1027        const int y=m_mpiInfo->rank/m_NX;
1028        for (int i1=-1; i1<2; i1++) {
1029            for (int i0=-1; i0<2; i0++) {
1030                // skip this rank
1031                if (i0==0 && i1==0)
1032                    continue;
1033                // location of neighbour rank
1034                const int nx=x+i0;
1035                const int ny=y+i1;
1036                if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {
1037                    neighbour.push_back(ny*m_NX+nx);
1038                    if (i0==0) {
1039                        // sharing top or bottom edge
1040                        const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1041                        const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);
1042                        offsetInShared.push_back(offsetInShared.back()+nDOF0);
1043                        for (dim_t i=0; i<nDOF0; i++, numShared++) {
1044                            sendShared.push_back(firstDOF+i);
1045                            recvShared.push_back(numDOF+numShared);
1046                            if (i>0)
1047                                colIndices[firstDOF+i-1].push_back(numShared);
1048                            colIndices[firstDOF+i].push_back(numShared);
1049                            if (i<nDOF0-1)
1050                                colIndices[firstDOF+i+1].push_back(numShared);
1051                            m_dofMap[firstNode+i]=numDOF+numShared;
1052                        }
1053                    } else if (i1==0) {
1054                        // sharing left or right edge
1055                        const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1056                        const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);
1057                        offsetInShared.push_back(offsetInShared.back()+nDOF1);
1058                        for (dim_t i=0; i<nDOF1; i++, numShared++) {
1059                            sendShared.push_back(firstDOF+i*nDOF0);
1060                            recvShared.push_back(numDOF+numShared);
1061                            if (i>0)
1062                                colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);
1063                            colIndices[firstDOF+i*nDOF0].push_back(numShared);
1064                            if (i<nDOF1-1)
1065                                colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1066                            m_dofMap[firstNode+i*m_N0]=numDOF+numShared;
1067                        }
1068                    } else {
1069                        // sharing a node
1070                        const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1071                        const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);
1072                        offsetInShared.push_back(offsetInShared.back()+1);
1073                        sendShared.push_back(dof);
1074                        recvShared.push_back(numDOF+numShared);
1075                        colIndices[dof].push_back(numShared);
1076                        m_dofMap[node]=numDOF+numShared;
1077                        ++numShared;
1078                    }
1079                }
1080          }          }
1081      }      }
     if (x>0) {  
         // left  
         index.push_back(node-1);  
         num++;  
     }  
1082    
1083      return num;      // create connector
1084  }      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1085                numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
1086                &offsetInShared[0], 1, 0, m_mpiInfo);
1087        Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
1088                numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
1089                &offsetInShared[0], 1, 0, m_mpiInfo);
1090        m_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1091        Paso_SharedComponents_free(snd_shcomp);
1092        Paso_SharedComponents_free(rcv_shcomp);
1093    
1094  //private      // create main and couple blocks
1095  void Rectangle::generateCouplePatterns(Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const      Paso_Pattern *mainPattern = createMainPattern();
1096  {      Paso_Pattern *colPattern, *rowPattern;
1097      IndexVector ptr(1,0);      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
     IndexVector index;  
     const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);  
     const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);  
     const IndexVector faces=getNumFacesPerBoundary();  
1098    
1099      // bottom edge      // allocate paso distribution
1100      dim_t n=0;      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
1101      if (faces[0] == 0) {              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
         index.push_back(2*(myN0+myN1+1));  
         index.push_back(2*(myN0+myN1+1)+1);  
         n+=2;  
         if (faces[2] == 0) {  
             index.push_back(0);  
             index.push_back(1);  
             index.push_back(2);  
             n+=3;  
         }  
     } else if (faces[2] == 0) {  
         index.push_back(1);  
         index.push_back(2);  
         n+=2;  
     }  
     // n=neighbours of bottom-left corner node  
     ptr.push_back(ptr.back()+n);  
     n=0;  
     if (faces[2] == 0) {  
         for (dim_t i=1; i<myN0-1; i++) {  
             index.push_back(i);  
             index.push_back(i+1);  
             index.push_back(i+2);  
             ptr.push_back(ptr.back()+3);  
         }  
         index.push_back(myN0-1);  
         index.push_back(myN0);  
         n+=2;  
         if (faces[1] == 0) {  
             index.push_back(myN0+1);  
             index.push_back(myN0+2);  
             index.push_back(myN0+3);  
             n+=3;  
         }  
     } else {  
         for (dim_t i=1; i<myN0-1; i++) {  
             ptr.push_back(ptr.back());  
         }  
         if (faces[1] == 0) {  
             index.push_back(myN0+2);  
             index.push_back(myN0+3);  
             n+=2;  
         }  
     }  
     // n=neighbours of bottom-right corner node  
     ptr.push_back(ptr.back()+n);  
1102    
1103      // 2nd row to 2nd last row      // finally create the system matrix
1104      for (dim_t i=1; i<myN1-1; i++) {      m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
1105          // left edge              distribution, distribution, mainPattern, colPattern, rowPattern,
1106          if (faces[0] == 0) {              m_connector, m_connector);
             index.push_back(2*(myN0+myN1+2)-i);  
             index.push_back(2*(myN0+myN1+2)-i-1);  
             index.push_back(2*(myN0+myN1+2)-i-2);  
             ptr.push_back(ptr.back()+3);  
         } else {  
             ptr.push_back(ptr.back());  
         }  
         for (dim_t j=1; j<myN0-1; j++) {  
             ptr.push_back(ptr.back());  
         }  
         // right edge  
         if (faces[1] == 0) {  
             index.push_back(myN0+i+1);  
             index.push_back(myN0+i+2);  
             index.push_back(myN0+i+3);  
             ptr.push_back(ptr.back()+3);  
         } else {  
             ptr.push_back(ptr.back());  
         }  
     }  
1107    
1108      // top edge      Paso_Distribution_free(distribution);
     n=0;  
     if (faces[0] == 0) {  
         index.push_back(2*myN0+myN1+5);  
         index.push_back(2*myN0+myN1+4);  
         n+=2;  
         if (faces[3] == 0) {  
             index.push_back(2*myN0+myN1+3);  
             index.push_back(2*myN0+myN1+2);  
             index.push_back(2*myN0+myN1+1);  
             n+=3;  
         }  
     } else if (faces[3] == 0) {  
         index.push_back(2*myN0+myN1+2);  
         index.push_back(2*myN0+myN1+1);  
         n+=2;  
     }  
     // n=neighbours of top-left corner node  
     ptr.push_back(ptr.back()+n);  
     n=0;  
     if (faces[3] == 0) {  
         for (dim_t i=1; i<myN0-1; i++) {  
             index.push_back(2*myN0+myN1+i+1);  
             index.push_back(2*myN0+myN1+i);  
             index.push_back(2*myN0+myN1+i-1);  
             ptr.push_back(ptr.back()+3);  
         }  
         index.push_back(myN0+myN1+4);  
         index.push_back(myN0+myN1+3);  
         n+=2;  
         if (faces[1] == 0) {  
             index.push_back(myN0+myN1+2);  
             index.push_back(myN0+myN1+1);  
             index.push_back(myN0+myN1);  
             n+=3;  
         }  
     } else {  
         for (dim_t i=1; i<myN0-1; i++) {  
             ptr.push_back(ptr.back());  
         }  
         if (faces[1] == 0) {  
             index.push_back(myN0+myN1+1);  
             index.push_back(myN0+myN1);  
             n+=2;  
         }  
     }  
     // n=neighbours of top-right corner node  
     ptr.push_back(ptr.back()+n);  
1109    
1110      dim_t M=ptr.size()-1;      // useful debug output
1111      map<index_t,index_t> idMap;      /*
1112      dim_t N=0;      cout << "--- rcv_shcomp ---" << endl;
1113      for (IndexVector::iterator it=index.begin(); it!=index.end(); it++) {      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
1114          if (idMap.find(*it)==idMap.end()) {      for (size_t i=0; i<neighbour.size(); i++) {
1115              idMap[*it]=N;          cout << "neighbor[" << i << "]=" << neighbour[i]
1116              *it=N++;              << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
1117          } else {      }
1118              *it=idMap[*it];      for (size_t i=0; i<recvShared.size(); i++) {
1119          }          cout << "shared[" << i << "]=" << recvShared[i] << endl;
1120      }      }
1121        cout << "--- snd_shcomp ---" << endl;
1122        for (size_t i=0; i<sendShared.size(); i++) {
1123            cout << "shared[" << i << "]=" << sendShared[i] << endl;
1124        }
1125        cout << "--- dofMap ---" << endl;
1126        for (size_t i=0; i<m_dofMap.size(); i++) {
1127            cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
1128        }
1129        cout << "--- colIndices ---" << endl;
1130        for (size_t i=0; i<colIndices.size(); i++) {
1131            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
1132        }
1133        */
1134    
1135      /*      /*
1136      cout << "--- colCouple_pattern ---" << endl;      cout << "--- main_pattern ---" << endl;
1137      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
1138      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<mainPattern->numOutput+1; i++) {
1139          cout << "ptr[" << i << "]=" << ptr[i] << endl;          cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
1140      }      }
1141      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
1142          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
1143      }      }
1144      */      */
1145    
1146      // now build the row couple pattern      /*
1147      IndexVector ptr2(1,0);      cout << "--- colCouple_pattern ---" << endl;
1148      IndexVector index2;      cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
1149      for (dim_t id=0; id<N; id++) {      for (size_t i=0; i<colPattern->numOutput+1; i++) {
1150          n=0;          cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
1151          for (dim_t i=0; i<M; i++) {      }
1152              for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {      for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
1153                  if (index[j]==id) {          cout << "index[" << i << "]=" << colPattern->index[i] << endl;
                     index2.push_back(i);  
                     n++;  
                     break;  
                 }  
             }  
         }  
         ptr2.push_back(ptr2.back()+n);  
1154      }      }
1155        */
1156    
1157      /*      /*
1158      cout << "--- rowCouple_pattern ---" << endl;      cout << "--- rowCouple_pattern ---" << endl;
1159      cout << "M=" << N << ", N=" << M << endl;      cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
1160      for (size_t i=0; i<ptr2.size(); i++) {      for (size_t i=0; i<rowPattern->numOutput+1; i++) {
1161          cout << "ptr[" << i << "]=" << ptr2[i] << endl;          cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
1162      }      }
1163      for (size_t i=0; i<index2.size(); i++) {      for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
1164          cout << "index[" << i << "]=" << index2[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
1165      }      }
1166      */      */
1167    
1168      // paso will manage the memory      Paso_Pattern_free(mainPattern);
1169      index_t* indexC = MEMALLOC(index.size(), index_t);      Paso_Pattern_free(colPattern);
1170      index_t* ptrC = MEMALLOC(ptr.size(), index_t);      Paso_Pattern_free(rowPattern);
     copy(index.begin(), index.end(), indexC);  
     copy(ptr.begin(), ptr.end(), ptrC);  
     *colPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);  
   
     // paso will manage the memory  
     indexC = MEMALLOC(index2.size(), index_t);  
     ptrC = MEMALLOC(ptr2.size(), index_t);  
     copy(index2.begin(), index2.end(), indexC);  
     copy(ptr2.begin(), ptr2.end(), ptrC);  
     *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);  
1171  }  }
1172    
1173  //protected  //protected
# Line 1512  void Rectangle::interpolateNodesOnFaces( Line 1332  void Rectangle::interpolateNodesOnFaces(
1332      }      }
1333  }  }
1334    
1335    //protected
1336    void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1337            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1338            const escript::Data& C, const escript::Data& D,
1339            const escript::Data& X, const escript::Data& Y,
1340            const escript::Data& d, const escript::Data& y) const
1341    {
1342        const double h0 = m_l0/m_gNE0;
1343        const double h1 = m_l1/m_gNE1;
1344        /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */
1345        const double w0 = -0.1555021169820365539*h1/h0;
1346        const double w1 = 0.041666666666666666667;
1347        const double w10 = -0.041666666666666666667*h0/h1;
1348        const double w11 = 0.1555021169820365539*h1/h0;
1349        const double w12 = 0.1555021169820365539*h0/h1;
1350        const double w13 = 0.01116454968463011277*h0/h1;
1351        const double w14 = 0.01116454968463011277*h1/h0;
1352        const double w15 = 0.041666666666666666667*h1/h0;
1353        const double w16 = -0.01116454968463011277*h0/h1;
1354        const double w17 = -0.1555021169820365539*h0/h1;
1355        const double w18 = -0.33333333333333333333*h1/h0;
1356        const double w19 = 0.25000000000000000000;
1357        const double w2 = -0.15550211698203655390;
1358        const double w20 = -0.25000000000000000000;
1359        const double w21 = 0.16666666666666666667*h0/h1;
1360        const double w22 = -0.16666666666666666667*h1/h0;
1361        const double w23 = -0.16666666666666666667*h0/h1;
1362        const double w24 = 0.33333333333333333333*h1/h0;
1363        const double w25 = 0.33333333333333333333*h0/h1;
1364        const double w26 = 0.16666666666666666667*h1/h0;
1365        const double w27 = -0.33333333333333333333*h0/h1;
1366        const double w28 = -0.032861463941450536761*h1;
1367        const double w29 = -0.032861463941450536761*h0;
1368        const double w3 = 0.041666666666666666667*h0/h1;
1369        const double w30 = -0.12264065304058601714*h1;
1370        const double w31 = -0.0023593469594139828636*h1;
1371        const double w32 = -0.008805202725216129906*h0;
1372        const double w33 = -0.008805202725216129906*h1;
1373        const double w34 = 0.032861463941450536761*h1;
1374        const double w35 = 0.008805202725216129906*h1;
1375        const double w36 = 0.008805202725216129906*h0;
1376        const double w37 = 0.0023593469594139828636*h1;
1377        const double w38 = 0.12264065304058601714*h1;
1378        const double w39 = 0.032861463941450536761*h0;
1379        const double w4 = 0.15550211698203655390;
1380        const double w40 = -0.12264065304058601714*h0;
1381        const double w41 = -0.0023593469594139828636*h0;
1382        const double w42 = 0.0023593469594139828636*h0;
1383        const double w43 = 0.12264065304058601714*h0;
1384        const double w44 = -0.16666666666666666667*h1;
1385        const double w45 = -0.083333333333333333333*h0;
1386        const double w46 = 0.083333333333333333333*h1;
1387        const double w47 = 0.16666666666666666667*h1;
1388        const double w48 = 0.083333333333333333333*h0;
1389        const double w49 = -0.16666666666666666667*h0;
1390        const double w5 = -0.041666666666666666667;
1391        const double w50 = 0.16666666666666666667*h0;
1392        const double w51 = -0.083333333333333333333*h1;
1393        const double w52 = 0.025917019497006092316*h0*h1;
1394        const double w53 = 0.0018607582807716854616*h0*h1;
1395        const double w54 = 0.0069444444444444444444*h0*h1;
1396        const double w55 = 0.09672363354357992482*h0*h1;
1397        const double w56 = 0.00049858867864229740201*h0*h1;
1398        const double w57 = 0.055555555555555555556*h0*h1;
1399        const double w58 = 0.027777777777777777778*h0*h1;
1400        const double w59 = 0.11111111111111111111*h0*h1;
1401        const double w6 = -0.01116454968463011277*h1/h0;
1402        const double w60 = -0.19716878364870322056*h1;
1403        const double w61 = -0.19716878364870322056*h0;
1404        const double w62 = -0.052831216351296779436*h0;
1405        const double w63 = -0.052831216351296779436*h1;
1406        const double w64 = 0.19716878364870322056*h1;
1407        const double w65 = 0.052831216351296779436*h1;
1408        const double w66 = 0.19716878364870322056*h0;
1409        const double w67 = 0.052831216351296779436*h0;
1410        const double w68 = -0.5*h1;
1411        const double w69 = -0.5*h0;
1412        const double w7 = 0.011164549684630112770;
1413        const double w70 = 0.5*h1;
1414        const double w71 = 0.5*h0;
1415        const double w72 = 0.1555021169820365539*h0*h1;
1416        const double w73 = 0.041666666666666666667*h0*h1;
1417        const double w74 = 0.01116454968463011277*h0*h1;
1418        const double w75 = 0.25*h0*h1;
1419        const double w8 = -0.011164549684630112770;
1420        const double w9 = -0.041666666666666666667*h1/h0;
1421        /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */
1422    
1423        rhs.requireWrite();
1424    #pragma omp parallel
1425        {
1426            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1427    #pragma omp for
1428                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
1429                    for (index_t k0=0; k0<m_NE0; ++k0)  {
1430                        bool add_EM_S=false;
1431                        bool add_EM_F=false;
1432                        vector<double> EM_S(4*4, 0);
1433                        vector<double> EM_F(4, 0);
1434                        const index_t e = k0 + m_NE0*k1;
1435                        /*** GENERATOR SNIP_PDE_SINGLE TOP */
1436                        ///////////////
1437                        // process A //
1438                        ///////////////
1439                        if (!A.isEmpty()) {
1440                            add_EM_S=true;
1441                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1442                            if (A.actsExpanded()) {
1443                                const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1444                                const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1445                                const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1446                                const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1447                                const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1448                                const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1449                                const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1450                                const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1451                                const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1452                                const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1453                                const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1454                                const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1455                                const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1456                                const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1457                                const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1458                                const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1459                                const register double tmp4_0 = A_10_1 + A_10_2;
1460                                const register double tmp12_0 = A_11_0 + A_11_2;
1461                                const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1462                                const register double tmp10_0 = A_01_3 + A_10_3;
1463                                const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
1464                                const register double tmp0_0 = A_01_0 + A_01_3;
1465                                const register double tmp13_0 = A_01_2 + A_10_1;
1466                                const register double tmp3_0 = A_00_2 + A_00_3;
1467                                const register double tmp11_0 = A_11_1 + A_11_3;
1468                                const register double tmp18_0 = A_01_1 + A_10_1;
1469                                const register double tmp1_0 = A_00_0 + A_00_1;
1470                                const register double tmp15_0 = A_01_1 + A_10_2;
1471                                const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1472                                const register double tmp16_0 = A_10_0 + A_10_3;
1473                                const register double tmp6_0 = A_01_3 + A_10_0;
1474                                const register double tmp17_0 = A_01_1 + A_01_2;
1475                                const register double tmp9_0 = A_01_0 + A_10_0;
1476                                const register double tmp7_0 = A_01_0 + A_10_3;
1477                                const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1478                                const register double tmp19_0 = A_01_2 + A_10_2;
1479                                const register double tmp14_1 = A_10_0*w8;
1480                                const register double tmp23_1 = tmp3_0*w14;
1481                                const register double tmp35_1 = A_01_0*w8;
1482                                const register double tmp54_1 = tmp13_0*w8;
1483                                const register double tmp20_1 = tmp9_0*w4;
1484                                const register double tmp25_1 = tmp12_0*w12;
1485                                const register double tmp2_1 = A_01_1*w4;
1486                                const register double tmp44_1 = tmp7_0*w7;
1487                                const register double tmp26_1 = tmp10_0*w4;
1488                                const register double tmp52_1 = tmp18_0*w8;
1489                                const register double tmp48_1 = A_10_1*w7;
1490                                const register double tmp46_1 = A_01_3*w8;
1491                                const register double tmp50_1 = A_01_0*w2;
1492                                const register double tmp8_1 = tmp4_0*w5;
1493                                const register double tmp56_1 = tmp19_0*w8;
1494                                const register double tmp9_1 = tmp2_0*w10;
1495                                const register double tmp19_1 = A_10_3*w2;
1496                                const register double tmp47_1 = A_10_2*w4;
1497                                const register double tmp16_1 = tmp3_0*w0;
1498                                const register double tmp18_1 = tmp1_0*w6;
1499                                const register double tmp31_1 = tmp11_0*w12;
1500                                const register double tmp55_1 = tmp15_0*w2;
1501                                const register double tmp39_1 = A_10_2*w7;
1502                                const register double tmp11_1 = tmp6_0*w7;
1503                                const register double tmp40_1 = tmp11_0*w17;
1504                                const register double tmp34_1 = tmp15_0*w8;
1505                                const register double tmp33_1 = tmp14_0*w5;
1506                                const register double tmp24_1 = tmp11_0*w13;
1507                                const register double tmp3_1 = tmp1_0*w0;
1508                                const register double tmp5_1 = tmp2_0*w3;
1509                                const register double tmp43_1 = tmp17_0*w5;
1510                                const register double tmp15_1 = A_01_2*w4;
1511                                const register double tmp53_1 = tmp19_0*w2;
1512                                const register double tmp27_1 = tmp3_0*w11;
1513                                const register double tmp32_1 = tmp13_0*w2;
1514                                const register double tmp10_1 = tmp5_0*w9;
1515                                const register double tmp37_1 = A_10_1*w4;
1516                                const register double tmp38_1 = tmp5_0*w15;
1517                                const register double tmp17_1 = A_01_1*w7;
1518                                const register double tmp12_1 = tmp7_0*w4;
1519                                const register double tmp22_1 = tmp10_0*w7;
1520                                const register double tmp57_1 = tmp18_0*w2;
1521                                const register double tmp28_1 = tmp9_0*w7;
1522                                const register double tmp29_1 = tmp1_0*w14;
1523                                const register double tmp51_1 = tmp11_0*w16;
1524                                const register double tmp42_1 = tmp12_0*w16;
1525                                const register double tmp49_1 = tmp12_0*w17;
1526                                const register double tmp21_1 = tmp1_0*w11;
1527                                const register double tmp1_1 = tmp0_0*w1;
1528                                const register double tmp45_1 = tmp6_0*w4;
1529                                const register double tmp7_1 = A_10_0*w2;
1530                                const register double tmp6_1 = tmp3_0*w6;
1531                                const register double tmp13_1 = tmp8_0*w1;
1532                                const register double tmp36_1 = tmp16_0*w1;
1533                                const register double tmp41_1 = A_01_3*w2;
1534                                const register double tmp30_1 = tmp12_0*w13;
1535                                const register double tmp4_1 = A_01_2*w7;
1536                                const register double tmp0_1 = A_10_3*w8;
1537                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
1538                                EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
1539                                EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
1540                                EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1541                                EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1542                                EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
1543                                EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
1544                                EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
1545                                EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1546                                EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
1547                                EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
1548                                EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
1549                                EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
1550                                EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
1551                                EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
1552                                EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
1553                            } else { /* constant data */
1554                                const register double A_00 = A_p[INDEX2(0,0,2)];
1555                                const register double A_01 = A_p[INDEX2(0,1,2)];
1556                                const register double A_10 = A_p[INDEX2(1,0,2)];
1557                                const register double A_11 = A_p[INDEX2(1,1,2)];
1558                                const register double tmp0_0 = A_01 + A_10;
1559                                const register double tmp0_1 = A_00*w18;
1560                                const register double tmp10_1 = A_01*w20;
1561                                const register double tmp12_1 = A_00*w26;
1562                                const register double tmp4_1 = A_00*w22;
1563                                const register double tmp8_1 = A_00*w24;
1564                                const register double tmp13_1 = A_10*w19;
1565                                const register double tmp9_1 = tmp0_0*w20;
1566                                const register double tmp3_1 = A_11*w21;
1567                                const register double tmp11_1 = A_11*w27;
1568                                const register double tmp1_1 = A_01*w19;
1569                                const register double tmp6_1 = A_11*w23;
1570                                const register double tmp7_1 = A_11*w25;
1571                                const register double tmp2_1 = A_10*w20;
1572                                const register double tmp5_1 = tmp0_0*w19;
1573                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1574                                EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1575                                EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1576                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1577                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1578                                EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1579                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1580                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1581                                EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1582                                EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1583                                EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1584                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1585                                EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1586                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1587                                EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1588                                EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1589                            }
1590                        }
1591                        ///////////////
1592                        // process B //
1593                        ///////////////
1594                        if (!B.isEmpty()) {
1595                            add_EM_S=true;
1596                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1597                            if (B.actsExpanded()) {
1598                                const register double B_0_0 = B_p[INDEX2(0,0,2)];
1599                                const register double B_1_0 = B_p[INDEX2(1,0,2)];
1600                                const register double B_0_1 = B_p[INDEX2(0,1,2)];
1601                                const register double B_1_1 = B_p[INDEX2(1,1,2)];
1602                                const register double B_0_2 = B_p[INDEX2(0,2,2)];
1603                                const register double B_1_2 = B_p[INDEX2(1,2,2)];
1604                                const register double B_0_3 = B_p[INDEX2(0,3,2)];
1605                                const register double B_1_3 = B_p[INDEX2(1,3,2)];
1606                                const register double tmp3_0 = B_0_0 + B_0_2;
1607                                const register double tmp1_0 = B_1_2 + B_1_3;
1608                                const register double tmp2_0 = B_0_1 + B_0_3;
1609                                const register double tmp0_0 = B_1_0 + B_1_1;
1610                                const register double tmp63_1 = B_1_1*w42;
1611                                const register double tmp79_1 = B_1_1*w40;
1612                                const register double tmp37_1 = tmp3_0*w35;
1613                                const register double tmp8_1 = tmp0_0*w32;
1614                                const register double tmp71_1 = B_0_1*w34;
1615                                const register double tmp19_1 = B_0_3*w31;
1616                                const register double tmp15_1 = B_0_3*w34;
1617                                const register double tmp9_1 = tmp3_0*w34;
1618                                const register double tmp35_1 = B_1_0*w36;
1619                                const register double tmp66_1 = B_0_3*w28;
1620                                const register double tmp28_1 = B_1_0*w42;
1621                                const register double tmp22_1 = B_1_0*w40;
1622                                const register double tmp16_1 = B_1_2*w29;
1623                                const register double tmp6_1 = tmp2_0*w35;
1624                                const register double tmp55_1 = B_1_3*w40;
1625                                const register double tmp50_1 = B_1_3*w42;
1626                                const register double tmp7_1 = tmp1_0*w29;
1627                                const register double tmp1_1 = tmp1_0*w32;
1628                                const register double tmp57_1 = B_0_3*w30;
1629                                const register double tmp18_1 = B_1_1*w32;
1630                                const register double tmp53_1 = B_1_0*w41;
1631                                const register double tmp61_1 = B_1_3*w36;
1632                                const register double tmp27_1 = B_0_3*w38;
1633                                const register double tmp64_1 = B_0_2*w30;
1634                                const register double tmp76_1 = B_0_1*w38;
1635                                const register double tmp39_1 = tmp2_0*w34;
1636                                const register double tmp62_1 = B_0_1*w31;
1637                                const register double tmp56_1 = B_0_0*w31;
1638                                const register double tmp49_1 = B_1_1*w36;
1639                                const register double tmp2_1 = B_0_2*w31;
1640                                const register double tmp23_1 = B_0_2*w33;
1641                                const register double tmp38_1 = B_1_1*w43;
1642                                const register double tmp74_1 = B_1_2*w41;
1643                                const register double tmp43_1 = B_1_1*w41;
1644                                const register double tmp58_1 = B_0_2*w28;
1645                                const register double tmp67_1 = B_0_0*w33;
1646                                const register double tmp33_1 = tmp0_0*w39;
1647                                const register double tmp4_1 = B_0_0*w28;
1648                                const register double tmp20_1 = B_0_0*w30;
1649                                const register double tmp13_1 = B_0_2*w38;
1650                                const register double tmp65_1 = B_1_2*w43;
1651                                const register double tmp0_1 = tmp0_0*w29;
1652                                const register double tmp41_1 = tmp3_0*w33;
1653                                const register double tmp73_1 = B_0_2*w37;
1654                                const register double tmp69_1 = B_0_0*w38;
1655                                const register double tmp48_1 = B_1_2*w39;
1656                                const register double tmp59_1 = B_0_1*w33;
1657                                const register double tmp17_1 = B_1_3*w41;
1658                                const register double tmp5_1 = B_0_3*w33;
1659                                const register double tmp3_1 = B_0_1*w30;
1660                                const register double tmp21_1 = B_0_1*w28;
1661                                const register double tmp42_1 = B_1_0*w29;
1662                                const register double tmp54_1 = B_1_2*w32;
1663                                const register double tmp60_1 = B_1_0*w39;
1664                                const register double tmp32_1 = tmp1_0*w36;
1665                                const register double tmp10_1 = B_0_1*w37;
1666                                const register double tmp14_1 = B_0_0*w35;
1667                                const register double tmp29_1 = B_0_1*w35;
1668                                const register double tmp26_1 = B_1_2*w36;
1669                                const register double tmp30_1 = B_1_3*w43;
1670                                const register double tmp70_1 = B_0_2*w35;
1671                                const register double tmp34_1 = B_1_3*w39;
1672                                const register double tmp51_1 = B_1_0*w43;
1673                                const register double tmp31_1 = B_0_2*w34;
1674                                const register double tmp45_1 = tmp3_0*w28;
1675                                const register double tmp11_1 = tmp1_0*w39;
1676                                const register double tmp52_1 = B_1_1*w29;
1677                                const register double tmp44_1 = B_1_3*w32;
1678                                const register double tmp25_1 = B_1_1*w39;
1679                                const register double tmp47_1 = tmp2_0*w33;
1680                                const register double tmp72_1 = B_1_3*w29;
1681                                const register double tmp40_1 = tmp2_0*w28;
1682                                const register double tmp46_1 = B_1_2*w40;
1683                                const register double tmp36_1 = B_1_2*w42;
1684                                const register double tmp24_1 = B_0_0*w37;
1685                                const register double tmp77_1 = B_0_3*w35;
1686                                const register double tmp68_1 = B_0_3*w37;
1687                                const register double tmp78_1 = B_0_0*w34;
1688                                const register double tmp12_1 = tmp0_0*w36;
1689                                const register double tmp75_1 = B_1_0*w32;
1690                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1691                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1692                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1693                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1694                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1695                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
1696                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1697                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1698                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1699                                EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1700                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1701                                EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1702                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1703                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1704                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
1705                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1706                            } else { /* constant data */
1707                                const register double B_0 = B_p[0];
1708                                const register double B_1 = B_p[1];
1709                                const register double tmp6_1 = B_1*w50;
1710                                const register double tmp1_1 = B_1*w45;
1711                                const register double tmp5_1 = B_1*w49;
1712                                const register double tmp4_1 = B_1*w48;
1713                                const register double tmp0_1 = B_0*w44;
1714                                const register double tmp2_1 = B_0*w46;
1715                                const register double tmp7_1 = B_0*w51;
1716                                const register double tmp3_1 = B_0*w47;
1717                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1718                                EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
1719                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1720                                EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1721                                EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
1722                                EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
1723                                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
1724                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
1725                                EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
1726                                EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
1727                                EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
1728                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
1729                                EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
1730                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1731                                EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
1732                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
1733                            }
1734                        }
1735                        ///////////////
1736                        // process C //
1737                        ///////////////
1738                        if (!C.isEmpty()) {
1739                            add_EM_S=true;
1740                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
1741                            if (C.actsExpanded()) {
1742                                const register double C_0_0 = C_p[INDEX2(0,0,2)];
1743                                const register double C_1_0 = C_p[INDEX2(1,0,2)];
1744                                const register double C_0_1 = C_p[INDEX2(0,1,2)];
1745                                const register double C_1_1 = C_p[INDEX2(1,1,2)];
1746                                const register double C_0_2 = C_p[INDEX2(0,2,2)];
1747                                const register double C_1_2 = C_p[INDEX2(1,2,2)];
1748                                const register double C_0_3 = C_p[INDEX2(0,3,2)];
1749                                const register double C_1_3 = C_p[INDEX2(1,3,2)];
1750                                const register double tmp2_0 = C_0_1 + C_0_3;
1751                                const register double tmp1_0 = C_1_2 + C_1_3;
1752                                const register double tmp3_0 = C_0_0 + C_0_2;
1753                                const register double tmp0_0 = C_1_0 + C_1_1;
1754                                const register double tmp64_1 = C_0_2*w30;
1755                                const register double tmp14_1 = C_0_2*w28;
1756                                const register double tmp19_1 = C_0_3*w31;
1757                                const register double tmp22_1 = C_1_0*w40;
1758                                const register double tmp37_1 = tmp3_0*w35;
1759                                const register double tmp29_1 = C_0_1*w35;
1760                                const register double tmp73_1 = C_0_2*w37;
1761                                const register double tmp74_1 = C_1_2*w41;
1762                                const register double tmp52_1 = C_1_3*w39;
1763                                const register double tmp25_1 = C_1_1*w39;
1764                                const register double tmp62_1 = C_0_1*w31;
1765                                const register double tmp79_1 = C_1_1*w40;
1766                                const register double tmp43_1 = C_1_1*w36;
1767                                const register double tmp27_1 = C_0_3*w38;
1768                                const register double tmp28_1 = C_1_0*w42;
1769                                const register double tmp63_1 = C_1_1*w42;
1770                                const register double tmp59_1 = C_0_3*w34;
1771                                const register double tmp72_1 = C_1_3*w29;
1772                                const register double tmp40_1 = tmp2_0*w35;
1773                                const register double tmp13_1 = C_0_3*w30;
1774                                const register double tmp51_1 = C_1_2*w40;
1775                                const register double tmp54_1 = C_1_2*w42;
1776                                const register double tmp12_1 = C_0_0*w31;
1777                                const register double tmp2_1 = tmp1_0*w32;
1778                                const register double tmp68_1 = C_0_2*w31;
1779                                const register double tmp75_1 = C_1_0*w32;
1780                                const register double tmp49_1 = C_1_1*w41;
1781                                const register double tmp4_1 = C_0_2*w35;
1782                                const register double tmp66_1 = C_0_3*w28;
1783                                const register double tmp56_1 = C_0_1*w37;
1784                                const register double tmp5_1 = C_0_1*w34;
1785                                const register double tmp38_1 = tmp2_0*w34;
1786                                const register double tmp76_1 = C_0_1*w38;
1787                                const register double tmp21_1 = C_0_1*w28;
1788                                const register double tmp69_1 = C_0_1*w30;
1789                                const register double tmp53_1 = C_1_0*w36;
1790                                const register double tmp42_1 = C_1_2*w39;
1791                                const register double tmp32_1 = tmp1_0*w29;
1792                                const register double tmp45_1 = C_1_0*w43;
1793                                const register double tmp33_1 = tmp0_0*w32;
1794                                const register double tmp35_1 = C_1_0*w41;
1795                                const register double tmp26_1 = C_1_2*w36;
1796                                const register double tmp67_1 = C_0_0*w33;
1797                                const register double tmp31_1 = C_0_2*w34;
1798                                const register double tmp20_1 = C_0_0*w30;
1799                                const register double tmp70_1 = C_0_0*w28;
1800                                const register double tmp8_1 = tmp0_0*w39;
1801                                const register double tmp30_1 = C_1_3*w43;
1802                                const register double tmp0_1 = tmp0_0*w29;
1803                                const register double tmp17_1 = C_1_3*w41;
1804                                const register double tmp58_1 = C_0_0*w35;
1805                                const register double tmp9_1 = tmp3_0*w33;
1806                                const register double tmp61_1 = C_1_3*w36;
1807                                const register double tmp41_1 = tmp3_0*w34;
1808                                const register double tmp50_1 = C_1_3*w32;
1809                                const register double tmp18_1 = C_1_1*w32;
1810                                const register double tmp6_1 = tmp1_0*w36;
1811                                const register double tmp3_1 = C_0_0*w38;
1812                                const register double tmp34_1 = C_1_1*w29;
1813                                const register double tmp77_1 = C_0_3*w35;
1814                                const register double tmp65_1 = C_1_2*w43;
1815                                const register double tmp71_1 = C_0_3*w33;
1816                                const register double tmp55_1 = C_1_1*w43;
1817                                const register double tmp46_1 = tmp3_0*w28;
1818                                const register double tmp24_1 = C_0_0*w37;
1819                                const register double tmp10_1 = tmp1_0*w39;
1820                                const register double tmp48_1 = C_1_0*w29;
1821                                const register double tmp15_1 = C_0_1*w33;
1822                                const register double tmp36_1 = C_1_2*w32;
1823                                const register double tmp60_1 = C_1_0*w39;
1824                                const register double tmp47_1 = tmp2_0*w33;
1825                                const register double tmp16_1 = C_1_2*w29;
1826                                const register double tmp1_1 = C_0_3*w37;
1827                                const register double tmp7_1 = tmp2_0*w28;
1828                                const register double tmp39_1 = C_1_3*w40;
1829                                const register double tmp44_1 = C_1_3*w42;
1830                                const register double tmp57_1 = C_0_2*w38;
1831                                const register double tmp78_1 = C_0_0*w34;
1832                                const register double tmp11_1 = tmp0_0*w36;
1833                                const register double tmp23_1 = C_0_2*w33;
1834                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1835                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1836                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1837                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1838                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1839                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
1840                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1841                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1842                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1843                                EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1844                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1845                                EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1846                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1847                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1848                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
1849                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1850                            } else { /* constant data */
1851                                const register double C_0 = C_p[0];
1852                                const register double C_1 = C_p[1];
1853                                const register double tmp1_1 = C_1*w45;
1854                                const register double tmp3_1 = C_0*w51;
1855                                const register double tmp4_1 = C_0*w44;
1856                                const register double tmp7_1 = C_0*w46;
1857                                const register double tmp5_1 = C_1*w49;
1858                                const register double tmp2_1 = C_1*w48;
1859                                const register double tmp0_1 = C_0*w47;
1860                                const register double tmp6_1 = C_1*w50;
1861                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1862                                EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
1863                                EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
1864                                EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1865                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
1866                                EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
1867                                EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
1868                                EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
1869                                EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
1870                                EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
1871                                EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
1872                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
1873                                EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
1874                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
1875                                EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
1876                                EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
1877                            }
1878                        }
1879                        ///////////////
1880                        // process D //
1881                        ///////////////
1882                        if (!D.isEmpty()) {
1883                            add_EM_S=true;
1884                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
1885                            if (D.actsExpanded()) {
1886                                const register double D_0 = D_p[0];
1887                                const register double D_1 = D_p[1];
1888                                const register double D_2 = D_p[2];
1889                                const register double D_3 = D_p[3];
1890                                const register double tmp4_0 = D_1 + D_3;
1891                                const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;
1892                                const register double tmp5_0 = D_0 + D_2;
1893                                const register double tmp0_0 = D_0 + D_1;
1894                                const register double tmp6_0 = D_0 + D_3;
1895                                const register double tmp1_0 = D_2 + D_3;
1896                                const register double tmp3_0 = D_1 + D_2;
1897                                const register double tmp16_1 = D_1*w56;
1898                                const register double tmp14_1 = tmp6_0*w54;
1899                                const register double tmp8_1 = D_3*w55;
1900                                const register double tmp2_1 = tmp2_0*w54;
1901                                const register double tmp12_1 = tmp5_0*w52;
1902                                const register double tmp4_1 = tmp0_0*w53;
1903                                const register double tmp3_1 = tmp1_0*w52;
1904                                const register double tmp13_1 = tmp4_0*w53;
1905                                const register double tmp10_1 = tmp4_0*w52;
1906                                const register double tmp1_1 = tmp1_0*w53;
1907                                const register double tmp7_1 = D_3*w56;
1908                                const register double tmp0_1 = tmp0_0*w52;
1909                                const register double tmp11_1 = tmp5_0*w53;
1910                                const register double tmp9_1 = D_0*w56;
1911                                const register double tmp5_1 = tmp3_0*w54;
1912                                const register double tmp18_1 = D_2*w56;
1913                                const register double tmp17_1 = D_1*w55;
1914                                const register double tmp6_1 = D_0*w55;
1915                                const register double tmp15_1 = D_2*w55;
1916                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1917                                EM_S[INDEX2(1,2,4)]+=tmp2_1;
1918                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1919                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
1920                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
1921                                EM_S[INDEX2(3,0,4)]+=tmp2_1;
1922                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
1923                                EM_S[INDEX2(2,1,4)]+=tmp2_1;
1924                                EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
1925                                EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
1926                                EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
1927                                EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
1928                                EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
1929                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
1930                                EM_S[INDEX2(0,3,4)]+=tmp2_1;
1931                                EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
1932                            } else { /* constant data */
1933                                const register double D_0 = D_p[0];
1934                                const register double tmp2_1 = D_0*w59;
1935                                const register double tmp1_1 = D_0*w58;
1936                                const register double tmp0_1 = D_0*w57;
1937                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
1938                                EM_S[INDEX2(1,2,4)]+=tmp1_1;
1939                                EM_S[INDEX2(3,2,4)]+=tmp0_1;
1940                                EM_S[INDEX2(0,0,4)]+=tmp2_1;
1941                                EM_S[INDEX2(3,3,4)]+=tmp2_1;
1942                                EM_S[INDEX2(3,0,4)]+=tmp1_1;
1943                                EM_S[INDEX2(3,1,4)]+=tmp0_1;
1944                                EM_S[INDEX2(2,1,4)]+=tmp1_1;
1945                                EM_S[INDEX2(0,2,4)]+=tmp0_1;
1946                                EM_S[INDEX2(2,0,4)]+=tmp0_1;
1947                                EM_S[INDEX2(1,3,4)]+=tmp0_1;
1948                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
1949                                EM_S[INDEX2(2,2,4)]+=tmp2_1;
1950                                EM_S[INDEX2(1,0,4)]+=tmp0_1;
1951                                EM_S[INDEX2(0,3,4)]+=tmp1_1;
1952                                EM_S[INDEX2(1,1,4)]+=tmp2_1;
1953                            }
1954                        }
1955                        ///////////////
1956                        // process X //
1957                        ///////////////
1958                        if (!X.isEmpty()) {
1959                            add_EM_F=true;
1960                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
1961                            if (X.actsExpanded()) {
1962                                const register double X_0_0 = X_p[INDEX2(0,0,2)];
1963                                const register double X_1_0 = X_p[INDEX2(1,0,2)];
1964                                const register double X_0_1 = X_p[INDEX2(0,1,2)];
1965                                const register double X_1_1 = X_p[INDEX2(1,1,2)];
1966                                const register double X_0_2 = X_p[INDEX2(0,2,2)];
1967                                const register double X_1_2 = X_p[INDEX2(1,2,2)];
1968                                const register double X_0_3 = X_p[INDEX2(0,3,2)];
1969                                const register double X_1_3 = X_p[INDEX2(1,3,2)];
1970                                const register double tmp1_0 = X_1_1 + X_1_3;
1971                                const register double tmp3_0 = X_0_0 + X_0_1;
1972                                const register double tmp2_0 = X_1_0 + X_1_2;
1973                                const register double tmp0_0 = X_0_2 + X_0_3;
1974                                const register double tmp8_1 = tmp2_0*w66;
1975                                const register double tmp5_1 = tmp3_0*w64;
1976                                const register double tmp14_1 = tmp0_0*w64;
1977                                const register double tmp3_1 = tmp3_0*w60;
1978                                const register double tmp9_1 = tmp3_0*w63;
1979                                const register double tmp13_1 = tmp3_0*w65;
1980                                const register double tmp12_1 = tmp1_0*w66;
1981                                const register double tmp10_1 = tmp0_0*w60;
1982                                const register double tmp2_1 = tmp2_0*w61;
1983                                const register double tmp6_1 = tmp2_0*w62;
1984                                const register double tmp4_1 = tmp0_0*w65;
1985                                const register double tmp11_1 = tmp1_0*w67;
1986                                const register double tmp1_1 = tmp1_0*w62;
1987                                const register double tmp7_1 = tmp1_0*w61;
1988                                const register double tmp0_1 = tmp0_0*w63;
1989                                const register double tmp15_1 = tmp2_0*w67;
1990                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1991                                EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
1992                                EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
1993                                EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1994                            } else { /* constant data */
1995                                const register double X_0 = X_p[0];
1996                                const register double X_1 = X_p[1];
1997                                const register double tmp3_1 = X_1*w71;
1998                                const register double tmp2_1 = X_0*w70;
1999                                const register double tmp1_1 = X_0*w68;
2000                                const register double tmp0_1 = X_1*w69;
2001                                EM_F[0]+=tmp0_1 + tmp1_1;
2002                                EM_F[1]+=tmp0_1 + tmp2_1;
2003                                EM_F[2]+=tmp1_1 + tmp3_1;
2004                                EM_F[3]+=tmp2_1 + tmp3_1;
2005                            }
2006                        }
2007                        ///////////////
2008                        // process Y //
2009                        ///////////////
2010                        if (!Y.isEmpty()) {
2011                            add_EM_F=true;
2012                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2013                            if (Y.actsExpanded()) {
2014                                const register double Y_0 = Y_p[0];
2015                                const register double Y_1 = Y_p[1];
2016                                const register double Y_2 = Y_p[2];
2017                                const register double Y_3 = Y_p[3];
2018                                const register double tmp0_0 = Y_1 + Y_2;
2019                                const register double tmp1_0 = Y_0 + Y_3;
2020                                const register double tmp9_1 = Y_0*w74;
2021                                const register double tmp4_1 = tmp1_0*w73;
2022                                const register double tmp5_1 = Y_2*w74;
2023                                const register double tmp7_1 = Y_1*w74;
2024                                const register double tmp6_1 = Y_2*w72;
2025                                const register double tmp2_1 = Y_3*w74;
2026                                const register double tmp8_1 = Y_3*w72;
2027                                const register double tmp3_1 = Y_1*w72;
2028                                const register double tmp0_1 = Y_0*w72;
2029                                const register double tmp1_1 = tmp0_0*w73;
2030                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2031                                EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2032                                EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2033                                EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2034                            } else { /* constant data */
2035                                const register double Y_0 = Y_p[0];
2036                                const register double tmp0_1 = Y_0*w75;
2037                                EM_F[0]+=tmp0_1;
2038                                EM_F[1]+=tmp0_1;
2039                                EM_F[2]+=tmp0_1;
2040                                EM_F[3]+=tmp0_1;
2041                            }
2042                        }
2043                        /* GENERATOR SNIP_PDE_SINGLE BOTTOM */
2044    
2045                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2046                        const index_t firstNode=m_N0*k1+k0;
2047                        IndexVector rowIndex;
2048                        rowIndex.push_back(m_dofMap[firstNode]);
2049                        rowIndex.push_back(m_dofMap[firstNode+1]);
2050                        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2051                        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2052                        if (add_EM_F) {
2053                            //cout << "-- AddtoRHS -- " << endl;
2054                            double *F_p=rhs.getSampleDataRW(0);
2055                            for (index_t i=0; i<rowIndex.size(); i++) {
2056                                if (rowIndex[i]<getNumDOF()) {
2057                                    F_p[rowIndex[i]]+=EM_F[i];
2058                                    //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;
2059                                }
2060                            }
2061                            //cout << "---"<<endl;
2062                        }
2063                        if (add_EM_S) {
2064                            //cout << "-- AddtoSystem -- " << endl;
2065                            addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2066                        }
2067    
2068                    } // end k0 loop
2069                } // end k1 loop
2070            } // end of coloring
2071        } // end of parallel region
2072    }
2073    
2074  } // end of namespace ripley  } // end of namespace ripley
2075    

Legend:
Removed from v.3733  
changed lines
  Added in v.3756

  ViewVC Help
Powered by ViewVC 1.1.26