/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

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

revision 3699 by caltinay, Thu Dec 1 22:59:42 2011 UTC revision 3759 by caltinay, Fri Jan 6 06:54:51 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 142  void Rectangle::dump(const string& fileN Line 162  void Rectangle::dump(const string& fileN
162      pair<double,double> ydy = getFirstCoordAndSpacing(1);      pair<double,double> ydy = getFirstCoordAndSpacing(1);
163  #pragma omp parallel  #pragma omp parallel
164      {      {
165  #pragma omp for  #pragma omp for nowait
166          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_N0; i0++) {
167              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=xdx.first+i0*xdx.second;
168          }          }
169  #pragma omp for  #pragma omp for nowait
170          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_N1; i1++) {
171              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=ydy.first+i1*ydy.second;
172          }          }
# 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 ReducedNodes: //FIXME: reduced
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:
251              return &m_elementId[0];              return &m_elementId[0];
252          case FaceElements:          case FaceElements:
253            case ReducedFaceElements:
254              return &m_faceId[0];              return &m_faceId[0];
255          default:          default:
256              break;              break;
# Line 232  const int* Rectangle::borrowSampleRefere Line 258  const int* Rectangle::borrowSampleRefere
258    
259      stringstream msg;      stringstream msg;
260      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs() not implemented for function space type "
261          << fsType;          << functionSpaceTypeAsString(fsType);
262      throw RipleyException(msg.str());      throw RipleyException(msg.str());
263  }  }
264    
265  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsType, index_t id) const
266  {  {
267  #ifdef ESYS_MPI      if (getMPISize()==1)
268      if (fsCode == Nodes) {          return true;
         const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];  
         const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;  
         return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);  
     } else  
         throw RipleyException("ownSample() only implemented for Nodes");  
 #else  
     return true;  
 #endif  
 }  
   
 void Rectangle::interpolateOnDomain(escript::Data& target,  
                                     const escript::Data& in) const  
 {  
     const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));  
     const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));  
     if (inDomain != *this)  
         throw RipleyException("Illegal domain of interpolant");  
     if (targetDomain != *this)  
         throw RipleyException("Illegal domain of interpolation target");  
269    
270      stringstream msg;      switch (fsType) {
     msg << "interpolateOnDomain() not implemented for function space "  
         << in.getFunctionSpace().getTypeCode() << " -> "  
         << target.getFunctionSpace().getTypeCode();  
   
     switch (in.getFunctionSpace().getTypeCode()) {  
271          case Nodes:          case Nodes:
272            case ReducedNodes: //FIXME: reduced
273                return (m_dofMap[id] < getNumDOF());
274          case DegreesOfFreedom:          case DegreesOfFreedom:
275              switch (target.getFunctionSpace().getTypeCode()) {          case ReducedDegreesOfFreedom:
276                  case DegreesOfFreedom:              return true;
277                      target=in;          case Elements:
278                      break;          case ReducedElements:
279                // check ownership of element's bottom left node
280                  case Elements:              return (m_dofMap[id%m_NE0+m_N0*(id/m_NE0)] < getNumDOF());
281                  {          case FaceElements:
282                      const double tmp0_2 = 0.62200846792814621559;          case ReducedFaceElements:
283                      const double tmp0_1 = 0.044658198738520451079;              {
284                      const double tmp0_0 = 0.16666666666666666667;                  // check ownership of face element's first node
285                      const dim_t numComp = in.getDataPointSize();                  const IndexVector faces = getNumFacesPerBoundary();
286                      escript::Data* inn=const_cast<escript::Data*>(&in);                  dim_t n=0;
287  #pragma omp parallel for                  for (size_t i=0; i<faces.size(); i++) {
288                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      n+=faces[i];
289                          for (index_t k0=0; k0 < m_NE0; ++k0) {                      if (id<n) {
290                              const register double* f_10 = inn->getSampleDataRO(INDEX2(k0+1,k1, m_N0));                          index_t k;
291                              const register double* f_11 = inn->getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                          if (i==1)
292                              const register double* f_01 = inn->getSampleDataRO(INDEX2(k0,k1+1, m_N0));                              k=m_N0-1;
293                              const register double* f_00 = inn->getSampleDataRO(INDEX2(k0,k1, m_N0));                          else if (i==3)
294                              double* o = target.getSampleDataRW(INDEX2(k0,k1,m_NE0));                              k=m_N0*(m_N1-1);
295                              for (index_t i=0; i < numComp; ++i) {                          else
296                                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);                              k=0;
297                                  o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);                          // determine whether to move right or up
298                                  o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);                          const index_t delta=(i/2==0 ? m_N0 : 1);
299                                  o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);                          return (m_dofMap[k+(id-n+faces[i])*delta] < getNumDOF());
300                              } // close component loop i                      }
                         } // close k0 loop  
                     } // close k1 loop  
                     break;  
301                  }                  }
302                    return false;
                 default:  
                     throw RipleyException(msg.str());  
303              }              }
             break;  
304          default:          default:
305              throw RipleyException(msg.str());              break;
306      }      }
307    
308        stringstream msg;
309        msg << "ownSample() not implemented for "
310            << functionSpaceTypeAsString(fsType);
311        throw RipleyException(msg.str());
312  }  }
313    
314  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const
                                                 bool reducedColOrder) const  
315  {  {
316      if (reducedRowOrder || reducedColOrder)      escript::Data& in = *const_cast<escript::Data*>(&cIn);
317          throw RipleyException("getPattern() not implemented for reduced order");      const dim_t numComp = in.getDataPointSize();
318        const double h0 = m_l0/m_gNE0;
319        const double h1 = m_l1/m_gNE1;
320        const double cx0 = -1./h0;
321        const double cx1 = -.78867513459481288225/h0;
322        const double cx2 = -.5/h0;
323        const double cx3 = -.21132486540518711775/h0;
324        const double cx4 = .21132486540518711775/h0;
325        const double cx5 = .5/h0;
326        const double cx6 = .78867513459481288225/h0;
327        const double cx7 = 1./h0;
328        const double cy0 = -1./h1;
329        const double cy1 = -.78867513459481288225/h1;
330        const double cy2 = -.5/h1;
331        const double cy3 = -.21132486540518711775/h1;
332        const double cy4 = .21132486540518711775/h1;
333        const double cy5 = .5/h1;
334        const double cy6 = .78867513459481288225/h1;
335        const double cy7 = 1./h1;
336    
337      // connector      if (out.getFunctionSpace().getTypeCode() == Elements) {
338      RankVector neighbour;          /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */
339      IndexVector offsetInShared(1,0);  #pragma omp parallel for
340      IndexVector sendShared, recvShared;          for (index_t k1=0; k1 < m_NE1; ++k1) {
341      const IndexVector faces=getNumFacesPerBoundary();              for (index_t k0=0; k0 < m_NE0; ++k0) {
342      const index_t left = (m_offset0==0 ? 0 : 1);                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
343      const index_t bottom = (m_offset1==0 ? 0 : 1);                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
344      // corner node from bottom-left                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
345      if (faces[0] == 0 && faces[2] == 0) {                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
346          neighbour.push_back(m_mpiInfo->rank-m_NX-1);                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
347          offsetInShared.push_back(offsetInShared.back()+1);                  for (index_t i=0; i < numComp; ++i) {
348          sendShared.push_back(m_nodeId[m_N0+left]);                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
349          recvShared.push_back(m_nodeId[0]);                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
350      }                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
351      // bottom edge                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
352      if (faces[2] == 0) {                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
353          neighbour.push_back(m_mpiInfo->rank-m_NX);                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
354          offsetInShared.push_back(offsetInShared.back()+m_N0-left);                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
355          for (dim_t i=left; i<m_N0; i++) {                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
356              // easy case, we know the neighbour id's                  } /* end of component loop i */
357              sendShared.push_back(m_nodeId[m_N0+i]);              } /* end of k0 loop */
358              recvShared.push_back(m_nodeId[i]);          } /* end of k1 loop */
359          }          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
360      }      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
361      // corner node from bottom-right          /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */
362      if (faces[1] == 0 && faces[2] == 0) {  #pragma omp parallel for
363          neighbour.push_back(m_mpiInfo->rank-m_NX+1);          for (index_t k1=0; k1 < m_NE1; ++k1) {
364          const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);              for (index_t k0=0; k0 < m_NE0; ++k0) {
365          const index_t N1=(neighbour.back()/m_NX == 0 ? m_N1 : m_N1-1);                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
366          const index_t first=m_nodeDistribution[neighbour.back()];                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
367          offsetInShared.push_back(offsetInShared.back()+1);                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
368          sendShared.push_back(m_nodeId[(bottom+1)*m_N0-1]);                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
369          recvShared.push_back(first+N0*(N1-1));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
370      }                  for (index_t i=0; i < numComp; ++i) {
371      // left edge                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
372      if (faces[0] == 0) {                      o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
373          neighbour.push_back(m_mpiInfo->rank-1);                  } /* end of component loop i */
374          offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);              } /* end of k0 loop */
375          for (dim_t i=bottom; i<m_N1; i++) {          } /* end of k1 loop */
376              // easy case, we know the neighbour id's          /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
377              sendShared.push_back(m_nodeId[i*m_N0+1]);      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
378              recvShared.push_back(m_nodeId[i*m_N0]);  #pragma omp parallel
379          }          {
380      }              /*** GENERATOR SNIP_GRAD_FACES TOP */
381      // right edge              if (m_faceOffset[0] > -1) {
382      if (faces[1] == 0) {  #pragma omp for nowait
383          neighbour.push_back(m_mpiInfo->rank+1);                  for (index_t k1=0; k1 < m_NE1; ++k1) {
384          const index_t rightN0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
385          const index_t first=m_nodeDistribution[neighbour.back()];                      const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
386          offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
387          for (dim_t i=bottom; i<m_N1; i++) {                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
388              sendShared.push_back(m_nodeId[(i+1)*m_N0-1]);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
389              recvShared.push_back(first+rightN0*(i-bottom));                      for (index_t i=0; i < numComp; ++i) {
390          }                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
391      }                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
392      // corner node from top-left                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
393      if (faces[0] == 0 && faces[3] == 0) {                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
394          neighbour.push_back(m_mpiInfo->rank+m_NX-1);                      } /* end of component loop i */
395          const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);                  } /* end of k1 loop */
396          const index_t first=m_nodeDistribution[neighbour.back()];              } /* end of face 0 */
397          offsetInShared.push_back(offsetInShared.back()+1);              if (m_faceOffset[1] > -1) {
398          sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+left]);  #pragma omp for nowait
399          recvShared.push_back(first+N0-1);                  for (index_t k1=0; k1 < m_NE1; ++k1) {
400      }                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
401      // top edge                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
402      if (faces[3] == 0) {                      const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
403          neighbour.push_back(m_mpiInfo->rank+m_NX);                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
404          const index_t first=m_nodeDistribution[neighbour.back()];                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
405          offsetInShared.push_back(offsetInShared.back()+m_N0-left);                      for (index_t i=0; i < numComp; ++i) {
406          for (dim_t i=left; i<m_N0; i++) {                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
407              sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+i]);                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
408              recvShared.push_back(first+i-left);                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
409          }                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
410      }                      } /* end of component loop i */
411      // corner node from top-right                  } /* end of k1 loop */
412      if (faces[1] == 0 && faces[3] == 0) {              } /* end of face 1 */
413          neighbour.push_back(m_mpiInfo->rank+m_NX+1);              if (m_faceOffset[2] > -1) {
414          const index_t first=m_nodeDistribution[neighbour.back()];  #pragma omp for nowait
415          offsetInShared.push_back(offsetInShared.back()+1);                  for (index_t k0=0; k0 < m_NE0; ++k0) {
416          sendShared.push_back(m_nodeId[m_N0*m_N1-1]);                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
417          recvShared.push_back(first);                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
418      }                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
419      const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
420      cout << "--- rcv_shcomp ---" << endl;                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
421      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;                      for (index_t i=0; i < numComp; ++i) {
422      for (size_t i=0; i<neighbour.size(); i++) {                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
423          cout << "neighbor[" << i << "]=" << neighbour[i]                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
424              << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
425      }                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
426      for (size_t i=0; i<recvShared.size(); i++) {                      } /* end of component loop i */
427          cout << "shared[" << i << "]=" << recvShared[i] << endl;                  } /* end of k0 loop */
428      }              } /* end of face 2 */
429      cout << "--- snd_shcomp ---" << endl;              if (m_faceOffset[3] > -1) {
430      for (size_t i=0; i<sendShared.size(); i++) {  #pragma omp for nowait
431          cout << "shared[" << i << "]=" << sendShared[i] << endl;                  for (index_t k0=0; k0 < m_NE0; ++k0) {
432                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
433                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
434                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
435                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
436                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
437                        for (index_t i=0; i < numComp; ++i) {
438                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
439                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
440                            o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
441                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
442                        } /* end of component loop i */
443                    } /* end of k0 loop */
444                } /* end of face 3 */
445                /* GENERATOR SNIP_GRAD_FACES BOTTOM */
446            } // end of parallel section
447        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
448    #pragma omp parallel
449            {
450                /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
451                if (m_faceOffset[0] > -1) {
452    #pragma omp for nowait
453                    for (index_t k1=0; k1 < m_NE1; ++k1) {
454                        const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
455                        const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
456                        const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
457                        const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
458                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
459                        for (index_t i=0; i < numComp; ++i) {
460                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
461                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
462                        } /* end of component loop i */
463                    } /* end of k1 loop */
464                } /* end of face 0 */
465                if (m_faceOffset[1] > -1) {
466    #pragma omp for nowait
467                    for (index_t k1=0; k1 < m_NE1; ++k1) {
468                        const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
469                        const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
470                        const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
471                        const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
472                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
473                        for (index_t i=0; i < numComp; ++i) {
474                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
475                            o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
476                        } /* end of component loop i */
477                    } /* end of k1 loop */
478                } /* end of face 1 */
479                if (m_faceOffset[2] > -1) {
480    #pragma omp for nowait
481                    for (index_t k0=0; k0 < m_NE0; ++k0) {
482                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
483                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
484                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
485                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
486                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
487                        for (index_t i=0; i < numComp; ++i) {
488                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
489                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
490                        } /* end of component loop i */
491                    } /* end of k0 loop */
492                } /* end of face 2 */
493                if (m_faceOffset[3] > -1) {
494    #pragma omp for nowait
495                    for (index_t k0=0; k0 < m_NE0; ++k0) {
496                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
497                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
498                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
499                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
500                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
501                        for (index_t i=0; i < numComp; ++i) {
502                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
503                            o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
504                        } /* end of component loop i */
505                    } /* end of k0 loop */
506                } /* end of face 3 */
507                /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */
508            } // end of parallel section
509        } else {
510            stringstream msg;
511            msg << "setToGradient() not implemented for "
512                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
513            throw RipleyException(msg.str());
514      }      }
515    }
516    
517      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(  void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
518              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],  {
519              &offsetInShared[0], 1, 0, m_mpiInfo);      escript::Data& in = *const_cast<escript::Data*>(&arg);
520      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      const dim_t numComp = in.getDataPointSize();
521              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],      const double h0 = m_l0/m_gNE0;
522              &offsetInShared[0], 1, 0, m_mpiInfo);      const double h1 = m_l1/m_gNE1;
523      Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);      if (arg.getFunctionSpace().getTypeCode() == Elements) {
524      Paso_SharedComponents_free(snd_shcomp);          const double w_0 = h0*h1/4.;
525      Paso_SharedComponents_free(rcv_shcomp);  #pragma omp parallel
526            {
527                vector<double> int_local(numComp, 0);
528    #pragma omp for nowait
529                for (index_t k1 = 0; k1 < m_NE1; ++k1) {
530                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
531                        const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
532                        for (index_t i=0; i < numComp; ++i) {
533                            const register double f_0 = f[INDEX2(i,0,numComp)];
534                            const register double f_1 = f[INDEX2(i,1,numComp)];
535                            const register double f_2 = f[INDEX2(i,2,numComp)];
536                            const register double f_3 = f[INDEX2(i,3,numComp)];
537                            int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
538                        }  /* end of component loop i */
539                    } /* end of k0 loop */
540                } /* end of k1 loop */
541    
542    #pragma omp critical
543                for (index_t i=0; i<numComp; i++)
544                    integrals[i]+=int_local[i];
545            } // end of parallel section
546        } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
547            const double w_0 = h0*h1;
548    #pragma omp parallel
549            {
550                vector<double> int_local(numComp, 0);
551    #pragma omp for nowait
552                for (index_t k1 = 0; k1 < m_NE1; ++k1) {
553                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
554                        const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
555                        for (index_t i=0; i < numComp; ++i) {
556                            int_local[i]+=f[i]*w_0;
557                        }  /* end of component loop i */
558                    } /* end of k0 loop */
559                } /* end of k1 loop */
560    
561    #pragma omp critical
562                for (index_t i=0; i<numComp; i++)
563                    integrals[i]+=int_local[i];
564            } // end of parallel section
565        } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
566            const double w_0 = h0/2.;
567            const double w_1 = h1/2.;
568    #pragma omp parallel
569            {
570                vector<double> int_local(numComp, 0);
571                if (m_faceOffset[0] > -1) {
572    #pragma omp for nowait
573                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
574                        const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
575                        for (index_t i=0; i < numComp; ++i) {
576                            const register double f_0 = f[INDEX2(i,0,numComp)];
577                            const register double f_1 = f[INDEX2(i,1,numComp)];
578                            int_local[i]+=(f_0+f_1)*w_1;
579                        }  /* end of component loop i */
580                    } /* end of k1 loop */
581                }
582    
583      // create patterns              if (m_faceOffset[1] > -1) {
584      dim_t M, N;  #pragma omp for nowait
585      IndexVector ptr(1,0);                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
586      IndexVector index;                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
587                        for (index_t i=0; i < numComp; ++i) {
588      // main pattern                          const register double f_0 = f[INDEX2(i,0,numComp)];
589      for (index_t i=0; i<numDOF; i++) {                          const register double f_1 = f[INDEX2(i,1,numComp)];
590          // always add the node itself                          int_local[i]+=(f_0+f_1)*w_1;
591          index.push_back(i);                      }  /* end of component loop i */
592          int num=insertNeighbours(index, i);                  } /* end of k1 loop */
593          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);  
594    
595      cout << "--- main_pattern ---" << endl;              if (m_faceOffset[2] > -1) {
596      cout << "M=" << M << ", N=" << N << endl;  #pragma omp for nowait
597      for (size_t i=0; i<ptr.size(); i++) {                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
598          cout << "ptr[" << i << "]=" << ptr[i] << endl;                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
599      }                      for (index_t i=0; i < numComp; ++i) {
600      for (size_t i=0; i<index.size(); i++) {                          const register double f_0 = f[INDEX2(i,0,numComp)];
601          cout << "index[" << i << "]=" << index[i] << endl;                          const register double f_1 = f[INDEX2(i,1,numComp)];
602                            int_local[i]+=(f_0+f_1)*w_0;
603                        }  /* end of component loop i */
604                    } /* end of k0 loop */
605                }
606    
607                if (m_faceOffset[3] > -1) {
608    #pragma omp for nowait
609                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
610                        const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
611                        for (index_t i=0; i < numComp; ++i) {
612                            const register double f_0 = f[INDEX2(i,0,numComp)];
613                            const register double f_1 = f[INDEX2(i,1,numComp)];
614                            int_local[i]+=(f_0+f_1)*w_0;
615                        }  /* end of component loop i */
616                    } /* end of k0 loop */
617                }
618    
619    #pragma omp critical
620                for (index_t i=0; i<numComp; i++)
621                    integrals[i]+=int_local[i];
622            } // end of parallel section
623        } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
624    #pragma omp parallel
625            {
626                vector<double> int_local(numComp, 0);
627                if (m_faceOffset[0] > -1) {
628    #pragma omp for nowait
629                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
630                        const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
631                        for (index_t i=0; i < numComp; ++i) {
632                            int_local[i]+=f[i]*h1;
633                        }  /* end of component loop i */
634                    } /* end of k1 loop */
635                }
636    
637                if (m_faceOffset[1] > -1) {
638    #pragma omp for nowait
639                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
640                        const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
641                        for (index_t i=0; i < numComp; ++i) {
642                            int_local[i]+=f[i]*h1;
643                        }  /* end of component loop i */
644                    } /* end of k1 loop */
645                }
646    
647                if (m_faceOffset[2] > -1) {
648    #pragma omp for nowait
649                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
650                        const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
651                        for (index_t i=0; i < numComp; ++i) {
652                            int_local[i]+=f[i]*h0;
653                        }  /* end of component loop i */
654                    } /* end of k0 loop */
655                }
656    
657                if (m_faceOffset[3] > -1) {
658    #pragma omp for nowait
659                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
660                        const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
661                        for (index_t i=0; i < numComp; ++i) {
662                            int_local[i]+=f[i]*h0;
663                        }  /* end of component loop i */
664                    } /* end of k0 loop */
665                }
666    
667    #pragma omp critical
668                for (index_t i=0; i<numComp; i++)
669                    integrals[i]+=int_local[i];
670            } // end of parallel section
671        } else {
672            stringstream msg;
673            msg << "setToIntegrals() not implemented for "
674                << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
675            throw RipleyException(msg.str());
676      }      }
677    }
678    
679      ptr.clear();  void Rectangle::setToNormal(escript::Data& out) const
680      index.clear();  {
681        if (out.getFunctionSpace().getTypeCode() == FaceElements) {
682    #pragma omp parallel
683            {
684                if (m_faceOffset[0] > -1) {
685    #pragma omp for nowait
686                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
687                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
688                        // set vector at two quadrature points
689                        *o++ = -1.;
690                        *o++ = 0.;
691                        *o++ = -1.;
692                        *o = 0.;
693                    }
694                }
695    
696      // column & row couple patterns              if (m_faceOffset[1] > -1) {
697      Paso_Pattern *colCouplePattern, *rowCouplePattern;  #pragma omp for nowait
698      generateCouplePatterns(&colCouplePattern, &rowCouplePattern);                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
699                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
700                        // set vector at two quadrature points
701                        *o++ = 1.;
702                        *o++ = 0.;
703                        *o++ = 1.;
704                        *o = 0.;
705                    }
706                }
707    
708      // allocate paso distribution              if (m_faceOffset[2] > -1) {
709      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,  #pragma omp for nowait
710              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
711                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
712                        // set vector at two quadrature points
713                        *o++ = 0.;
714                        *o++ = -1.;
715                        *o++ = 0.;
716                        *o = -1.;
717                    }
718                }
719    
720      Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(              if (m_faceOffset[3] > -1) {
721              MATRIX_FORMAT_DEFAULT, distribution, distribution,  #pragma omp for nowait
722              mainPattern, colCouplePattern, rowCouplePattern,                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
723              connector, connector);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
724      Paso_Pattern_free(mainPattern);                      // set vector at two quadrature points
725      Paso_Pattern_free(colCouplePattern);                      *o++ = 0.;
726      Paso_Pattern_free(rowCouplePattern);                      *o++ = 1.;
727      Paso_Distribution_free(distribution);                      *o++ = 0.;
728      return pattern;                      *o = 1.;
729                    }
730                }
731            } // end of parallel section
732        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
733    #pragma omp parallel
734            {
735                if (m_faceOffset[0] > -1) {
736    #pragma omp for nowait
737                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
738                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
739                        *o++ = -1.;
740                        *o = 0.;
741                    }
742                }
743    
744                if (m_faceOffset[1] > -1) {
745    #pragma omp for nowait
746                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
747                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
748                        *o++ = 1.;
749                        *o = 0.;
750                    }
751                }
752    
753                if (m_faceOffset[2] > -1) {
754    #pragma omp for nowait
755                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
756                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
757                        *o++ = 0.;
758                        *o = -1.;
759                    }
760                }
761    
762                if (m_faceOffset[3] > -1) {
763    #pragma omp for nowait
764                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
765                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
766                        *o++ = 0.;
767                        *o = 1.;
768                    }
769                }
770            } // end of parallel section
771    
772        } else {
773            stringstream msg;
774            msg << "setToNormal() not implemented for "
775                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
776            throw RipleyException(msg.str());
777        }
778    }
779    
780    Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
781                                                    bool reducedColOrder) const
782    {
783        if (reducedRowOrder || reducedColOrder)
784            throw RipleyException("getPattern() not implemented for reduced order");
785    
786        return m_pattern;
787  }  }
788    
789  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
# Line 535  pair<double,double> Rectangle::getFirstC Line 848  pair<double,double> Rectangle::getFirstC
848  }  }
849    
850  //protected  //protected
851    dim_t Rectangle::getNumDOF() const
852    {
853        return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;
854    }
855    
856    //protected
857  dim_t Rectangle::getNumFaceElements() const  dim_t Rectangle::getNumFaceElements() const
858  {  {
859      const IndexVector faces = getNumFacesPerBoundary();      const IndexVector faces = getNumFacesPerBoundary();
# Line 567  void Rectangle::assembleCoordinates(escr Line 886  void Rectangle::assembleCoordinates(escr
886      }      }
887  }  }
888    
889    //protected
890    dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
891    {
892        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
893        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
894        const int x=node%nDOF0;
895        const int y=node/nDOF0;
896        dim_t num=0;
897        // loop through potential neighbours and add to index if positions are
898        // within bounds
899        for (int i1=-1; i1<2; i1++) {
900            for (int i0=-1; i0<2; i0++) {
901                // skip node itself
902                if (i0==0 && i1==0)
903                    continue;
904                // location of neighbour node
905                const int nx=x+i0;
906                const int ny=y+i1;
907                if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {
908                    index.push_back(ny*nDOF0+nx);
909                    num++;
910                }
911            }
912        }
913    
914        return num;
915    }
916    
917    //protected
918    void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const
919    {
920        const dim_t numComp = in.getDataPointSize();
921        out.requireWrite();
922    
923        const index_t left = (m_offset0==0 ? 0 : 1);
924        const index_t bottom = (m_offset1==0 ? 0 : 1);
925        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
926        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
927    #pragma omp parallel for
928        for (index_t i=0; i<nDOF1; i++) {
929            for (index_t j=0; j<nDOF0; j++) {
930                const index_t n=j+left+(i+bottom)*m_N0;
931                const double* src=in.getSampleDataRO(n);
932                copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
933            }
934        }
935    }
936    
937    //protected
938    void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const
939    {
940        const dim_t numComp = in.getDataPointSize();
941        Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
942        in.requireWrite();
943        Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));
944    
945        const dim_t numDOF = getNumDOF();
946        out.requireWrite();
947        const double* buffer = Paso_Coupler_finishCollect(coupler);
948    
949    #pragma omp parallel for
950        for (index_t i=0; i<getNumNodes(); i++) {
951            const double* src=(m_dofMap[i]<numDOF ?
952                    in.getSampleDataRO(m_dofMap[i])
953                    : &buffer[(m_dofMap[i]-numDOF)*numComp]);
954            copy(src, src+numComp, out.getSampleDataRW(i));
955        }
956    }
957    
958  //private  //private
959  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
960  {  {
961      // identifiers are ordered from left to right, bottom to top on each rank,      // identifiers are ordered from left to right, bottom to top globablly.
     // except for the shared nodes which are owned by the rank below / to the  
     // left of the current rank  
962    
963      // build node distribution vector first.      // build node distribution vector first.
     // m_nodeDistribution[i] is the first node id on rank i, that is  
964      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
965      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
966      m_nodeDistribution[1]=getNumNodes();      const dim_t numDOF=getNumDOF();
967      for (dim_t k=1; k<m_mpiInfo->size-1; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
968          const index_t x=k%m_NX;          m_nodeDistribution[k]=k*numDOF;
         const index_t y=k/m_NX;  
         index_t numNodes=getNumNodes();  
         if (x>0)  
             numNodes-=m_N1;  
         if (y>0)  
             numNodes-=m_N0;  
         if (x>0 && y>0)  
             numNodes++; // subtracted corner twice -> fix that  
         m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;  
969      }      }
970      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
   
971      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
972        m_dofId.resize(numDOF);
973        m_elementId.resize(getNumElements());
974        m_faceId.resize(getNumFaceElements());
975    
976    #pragma omp parallel
977        {
978            // nodes
979    #pragma omp for nowait
980            for (dim_t i1=0; i1<m_N1; i1++) {
981                for (dim_t i0=0; i0<m_N0; i0++) {
982                    m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;
983                }
984            }
985    
986            // degrees of freedom
987    #pragma omp for nowait
988            for (dim_t k=0; k<numDOF; k++)
989                m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;
990    
991            // elements
992    #pragma omp for nowait
993            for (dim_t i1=0; i1<m_NE1; i1++) {
994                for (dim_t i0=0; i0<m_NE0; i0++) {
995                    m_elementId[i0+i1*m_NE0]=(m_offset1+i1)*m_gNE0+m_offset0+i0;
996                }
997            }
998    
999      // the bottom row and left column are not owned by this rank so the          // face elements
1000      // identifiers need to be computed accordingly  #pragma omp for
1001            for (dim_t k=0; k<getNumFaceElements(); k++)
1002                m_faceId[k]=k;
1003        } // end parallel section
1004    
1005        m_nodeTags.assign(getNumNodes(), 0);
1006        updateTagsInUse(Nodes);
1007    
1008        m_elementTags.assign(getNumElements(), 0);
1009        updateTagsInUse(Elements);
1010    
1011        // generate face offset vector and set face tags
1012        const IndexVector facesPerEdge = getNumFacesPerBoundary();
1013        const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
1014        const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
1015        m_faceOffset.assign(facesPerEdge.size(), -1);
1016        m_faceTags.clear();
1017        index_t offset=0;
1018        for (size_t i=0; i<facesPerEdge.size(); i++) {
1019            if (facesPerEdge[i]>0) {
1020                m_faceOffset[i]=offset;
1021                offset+=facesPerEdge[i];
1022                m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);
1023            }
1024        }
1025        setTagMap("left", LEFT);
1026        setTagMap("right", RIGHT);
1027        setTagMap("bottom", BOTTOM);
1028        setTagMap("top", TOP);
1029        updateTagsInUse(FaceElements);
1030    }
1031    
1032    //private
1033    void Rectangle::createPattern()
1034    {
1035        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1036        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1037      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1038      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1039      if (left>0) {  
1040          const int neighbour=m_mpiInfo->rank-1;      // populate node->DOF mapping with own degrees of freedom.
1041          const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      // The rest is assigned in the loop further down
1042        m_dofMap.assign(getNumNodes(), 0);
1043  #pragma omp parallel for  #pragma omp parallel for
1044          for (dim_t i1=bottom; i1<m_N1; i1++) {      for (index_t i=bottom; i<m_N1; i++) {
1045              m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]          for (index_t j=left; j<m_N0; j++) {
1046                  + (i1-bottom+1)*leftN0-1;              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
1047          }          }
1048      }      }
1049      if (bottom>0) {  
1050          const int neighbour=m_mpiInfo->rank-m_NX;      // build list of shared components and neighbours by looping through
1051          const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      // all potential neighbouring ranks and checking if positions are
1052          const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);      // within bounds
1053  #pragma omp parallel for      const dim_t numDOF=nDOF0*nDOF1;
1054          for (dim_t i0=left; i0<m_N0; i0++) {      vector<IndexVector> colIndices(numDOF); // for the couple blocks
1055              m_nodeId[i0]=m_nodeDistribution[neighbour]      RankVector neighbour;
1056                  + (bottomN1-1)*bottomN0 + i0 - left;      IndexVector offsetInShared(1,0);
1057        IndexVector sendShared, recvShared;
1058        int numShared=0;
1059        const int x=m_mpiInfo->rank%m_NX;
1060        const int y=m_mpiInfo->rank/m_NX;
1061        for (int i1=-1; i1<2; i1++) {
1062            for (int i0=-1; i0<2; i0++) {
1063                // skip this rank
1064                if (i0==0 && i1==0)
1065                    continue;
1066                // location of neighbour rank
1067                const int nx=x+i0;
1068                const int ny=y+i1;
1069                if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {
1070                    neighbour.push_back(ny*m_NX+nx);
1071                    if (i0==0) {
1072                        // sharing top or bottom edge
1073                        const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1074                        const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);
1075                        offsetInShared.push_back(offsetInShared.back()+nDOF0);
1076                        for (dim_t i=0; i<nDOF0; i++, numShared++) {
1077                            sendShared.push_back(firstDOF+i);
1078                            recvShared.push_back(numDOF+numShared);
1079                            if (i>0)
1080                                colIndices[firstDOF+i-1].push_back(numShared);
1081                            colIndices[firstDOF+i].push_back(numShared);
1082                            if (i<nDOF0-1)
1083                                colIndices[firstDOF+i+1].push_back(numShared);
1084                            m_dofMap[firstNode+i]=numDOF+numShared;
1085                        }
1086                    } else if (i1==0) {
1087                        // sharing left or right edge
1088                        const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1089                        const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);
1090                        offsetInShared.push_back(offsetInShared.back()+nDOF1);
1091                        for (dim_t i=0; i<nDOF1; i++, numShared++) {
1092                            sendShared.push_back(firstDOF+i*nDOF0);
1093                            recvShared.push_back(numDOF+numShared);
1094                            if (i>0)
1095                                colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);
1096                            colIndices[firstDOF+i*nDOF0].push_back(numShared);
1097                            if (i<nDOF1-1)
1098                                colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1099                            m_dofMap[firstNode+i*m_N0]=numDOF+numShared;
1100                        }
1101                    } else {
1102                        // sharing a node
1103                        const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1104                        const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);
1105                        offsetInShared.push_back(offsetInShared.back()+1);
1106                        sendShared.push_back(dof);
1107                        recvShared.push_back(numDOF+numShared);
1108                        colIndices[dof].push_back(numShared);
1109                        m_dofMap[node]=numDOF+numShared;
1110                        ++numShared;
1111                    }
1112                }
1113          }          }
1114      }      }
1115      if (left>0 && bottom>0) {  
1116          const int neighbour=m_mpiInfo->rank-m_NX-1;      // create connector
1117          const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1118          const index_t N1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
1119          m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1-1;              &offsetInShared[0], 1, 0, m_mpiInfo);
1120        Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
1121                numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
1122                &offsetInShared[0], 1, 0, m_mpiInfo);
1123        m_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1124        Paso_SharedComponents_free(snd_shcomp);
1125        Paso_SharedComponents_free(rcv_shcomp);
1126    
1127        // create main and couple blocks
1128        Paso_Pattern *mainPattern = createMainPattern();
1129        Paso_Pattern *colPattern, *rowPattern;
1130        createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
1131    
1132        // allocate paso distribution
1133        Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
1134                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
1135    
1136        // finally create the system matrix
1137        m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
1138                distribution, distribution, mainPattern, colPattern, rowPattern,
1139                m_connector, m_connector);
1140    
1141        Paso_Distribution_free(distribution);
1142    
1143        // useful debug output
1144        /*
1145        cout << "--- rcv_shcomp ---" << endl;
1146        cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
1147        for (size_t i=0; i<neighbour.size(); i++) {
1148            cout << "neighbor[" << i << "]=" << neighbour[i]
1149                << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
1150        }
1151        for (size_t i=0; i<recvShared.size(); i++) {
1152            cout << "shared[" << i << "]=" << recvShared[i] << endl;
1153      }      }
1154        cout << "--- snd_shcomp ---" << endl;
1155        for (size_t i=0; i<sendShared.size(); i++) {
1156            cout << "shared[" << i << "]=" << sendShared[i] << endl;
1157        }
1158        cout << "--- dofMap ---" << endl;
1159        for (size_t i=0; i<m_dofMap.size(); i++) {
1160            cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
1161        }
1162        cout << "--- colIndices ---" << endl;
1163        for (size_t i=0; i<colIndices.size(); i++) {
1164            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
1165        }
1166        */
1167    
1168      // the rest of the id's are contiguous      /*
1169      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];      cout << "--- main_pattern ---" << endl;
1170  #pragma omp parallel for      cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
1171      for (dim_t i1=bottom; i1<m_N1; i1++) {      for (size_t i=0; i<mainPattern->numOutput+1; i++) {
1172          for (dim_t i0=left; i0<m_N0; i0++) {          cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
1173              m_nodeId[i0+i1*m_N0] = firstId+i0-left+(i1-bottom)*(m_N0-left);      }
1174          }      for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
1175            cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
1176      }      }
1177        */
1178    
1179      // elements      /*
1180      m_elementId.resize(getNumElements());      cout << "--- colCouple_pattern ---" << endl;
1181  #pragma omp parallel for      cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
1182      for (dim_t k=0; k<getNumElements(); k++) {      for (size_t i=0; i<colPattern->numOutput+1; i++) {
1183          m_elementId[k]=k;          cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
1184        }
1185        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
1186            cout << "index[" << i << "]=" << colPattern->index[i] << endl;
1187      }      }
1188        */
1189    
1190      // face elements      /*
1191      m_faceId.resize(getNumFaceElements());      cout << "--- rowCouple_pattern ---" << endl;
1192  #pragma omp parallel for      cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
1193      for (dim_t k=0; k<getNumFaceElements(); k++) {      for (size_t i=0; i<rowPattern->numOutput+1; i++) {
1194          m_faceId[k]=k;          cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
1195        }
1196        for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
1197            cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
1198      }      }
1199        */
1200    
1201        Paso_Pattern_free(mainPattern);
1202        Paso_Pattern_free(colPattern);
1203        Paso_Pattern_free(rowPattern);
1204  }  }
1205    
1206  //private  //protected
1207  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const  void Rectangle::interpolateNodesOnElements(escript::Data& out,
1208                                            escript::Data& in, bool reduced) const
1209  {  {
1210      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);      const dim_t numComp = in.getDataPointSize();
1211      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);      if (reduced) {
1212      const int x=node%myN0;          /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */
1213      const int y=node/myN0;          const double c0 = .25;
1214      int num=0;  #pragma omp parallel for
1215      if (y>0) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1216          if (x>0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
1217              // bottom-left                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1218              index.push_back(node-myN0-1);                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1219              num++;                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1220          }                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1221          // bottom                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1222          index.push_back(node-myN0);                  for (index_t i=0; i < numComp; ++i) {
1223          num++;                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1224          if (x<myN0-1) {                  } /* end of component loop i */
1225              // bottom-right              } /* end of k0 loop */
1226              index.push_back(node-myN0+1);          } /* end of k1 loop */
1227              num++;          /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */
1228          }      } else {
1229      }          /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
1230      if (x<myN0-1) {          const double c0 = .16666666666666666667;
1231          // right          const double c1 = .044658198738520451079;
1232          index.push_back(node+1);          const double c2 = .62200846792814621559;
1233          num++;  #pragma omp parallel for
1234          if (y<myN1-1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1235              // top-right              for (index_t k0=0; k0 < m_NE0; ++k0) {
1236              index.push_back(node+myN0+1);                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1237              num++;                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1238          }                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1239      }                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1240      if (y<myN1-1) {                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1241          // top                  for (index_t i=0; i < numComp; ++i) {
1242          index.push_back(node+myN0);                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);
1243          num++;                      o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);
1244          if (x>0) {                      o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);
1245              // top-left                      o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);
1246              index.push_back(node+myN0-1);                  } /* end of component loop i */
1247              num++;              } /* end of k0 loop */
1248          }          } /* end of k1 loop */
1249      }          /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
     if (x>0) {  
         // left  
         index.push_back(node-1);  
         num++;  
1250      }      }
   
     return num;  
1251  }  }
1252    
1253  //private  //protected
1254  void Rectangle::generateCouplePatterns(Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const  void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,
1255                                            bool reduced) const
1256  {  {
1257      IndexVector ptr(1,0);      const dim_t numComp = in.getDataPointSize();
1258      IndexVector index;      if (reduced) {
1259      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);          const double c0 = .5;
1260      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);  #pragma omp parallel
1261      const IndexVector faces=getNumFacesPerBoundary();          {
1262                /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */
1263      // bottom edge              if (m_faceOffset[0] > -1) {
1264      dim_t n=0;  #pragma omp for nowait
1265      if (faces[0] == 0) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1266          index.push_back(2*(myN0+myN1+1));                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1267          index.push_back(2*(myN0+myN1+1)+1);                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1268          n+=2;                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1269          if (faces[2] == 0) {                      for (index_t i=0; i < numComp; ++i) {
1270              index.push_back(0);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
1271              index.push_back(1);                      } /* end of component loop i */
1272              index.push_back(2);                  } /* end of k1 loop */
1273              n+=3;              } /* end of face 0 */
1274          }              if (m_faceOffset[1] > -1) {
1275      } else if (faces[2] == 0) {  #pragma omp for nowait
1276          index.push_back(1);                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1277          index.push_back(2);                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1278          n+=2;                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1279      }                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1280      // n=neighbours of bottom-left corner node                      for (index_t i=0; i < numComp; ++i) {
1281      ptr.push_back(ptr.back()+n);                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
1282      n=0;                      } /* end of component loop i */
1283      if (faces[2] == 0) {                  } /* end of k1 loop */
1284          for (dim_t i=1; i<myN0-1; i++) {              } /* end of face 1 */
1285              index.push_back(i);              if (m_faceOffset[2] > -1) {
1286              index.push_back(i+1);  #pragma omp for nowait
1287              index.push_back(i+2);                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1288              ptr.push_back(ptr.back()+3);                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1289          }                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1290          index.push_back(myN0-1);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1291          index.push_back(myN0);                      for (index_t i=0; i < numComp; ++i) {
1292          n+=2;                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
1293          if (faces[1] == 0) {                      } /* end of component loop i */
1294              index.push_back(myN0+1);                  } /* end of k0 loop */
1295              index.push_back(myN0+2);              } /* end of face 2 */
1296              index.push_back(myN0+3);              if (m_faceOffset[3] > -1) {
1297              n+=3;  #pragma omp for nowait
1298          }                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1299      } else {                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1300          for (dim_t i=1; i<myN0-1; i++) {                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1301              ptr.push_back(ptr.back());                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1302          }                      for (index_t i=0; i < numComp; ++i) {
1303          if (faces[1] == 0) {                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1304              index.push_back(myN0+2);                      } /* end of component loop i */
1305              index.push_back(myN0+3);                  } /* end of k0 loop */
1306              n+=2;              } /* end of face 3 */
1307          }              /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */
1308      }          } // end of parallel section
     // n=neighbours of bottom-right corner node  
     ptr.push_back(ptr.back()+n);  
   
     // 2nd row to 2nd last row  
     for (dim_t i=1; i<myN1-1; i++) {  
         // left edge  
         if (faces[0] == 0) {  
             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());  
         }  
     }  
   
     // top edge  
     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;  
         }  
1309      } else {      } else {
1310          for (dim_t i=1; i<myN0-1; i++) {          const double c0 = 0.21132486540518711775;
1311              ptr.push_back(ptr.back());          const double c1 = 0.78867513459481288225;
1312          }  #pragma omp parallel
1313          if (faces[1] == 0) {          {
1314              index.push_back(myN0+myN1+1);              /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */
1315              index.push_back(myN0+myN1);              if (m_faceOffset[0] > -1) {
1316              n+=2;  #pragma omp for nowait
1317          }                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1318      }                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1319      // n=neighbours of top-right corner node                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1320      ptr.push_back(ptr.back()+n);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1321                        for (index_t i=0; i < numComp; ++i) {
1322      dim_t M=ptr.size()-1;                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;
1323      map<index_t,index_t> idMap;                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_01[i]*c1;
1324      dim_t N=0;                      } /* end of component loop i */
1325      for (IndexVector::iterator it=index.begin(); it!=index.end(); it++) {                  } /* end of k1 loop */
1326          if (idMap.find(*it)==idMap.end()) {              } /* end of face 0 */
1327              idMap[*it]=N;              if (m_faceOffset[1] > -1) {
1328              *it=N++;  #pragma omp for nowait
1329          } else {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1330              *it=idMap[*it];                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1331          }                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1332                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1333                        for (index_t i=0; i < numComp; ++i) {
1334                            o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;
1335                            o[INDEX2(i,numComp,1)] = f_10[i]*c0 + f_11[i]*c1;
1336                        } /* end of component loop i */
1337                    } /* end of k1 loop */
1338                } /* end of face 1 */
1339                if (m_faceOffset[2] > -1) {
1340    #pragma omp for nowait
1341                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1342                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1343                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1344                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1345                        for (index_t i=0; i < numComp; ++i) {
1346                            o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;
1347                            o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_10[i]*c1;
1348                        } /* end of component loop i */
1349                    } /* end of k0 loop */
1350                } /* end of face 2 */
1351                if (m_faceOffset[3] > -1) {
1352    #pragma omp for nowait
1353                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1354                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1355                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1356                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1357                        for (index_t i=0; i < numComp; ++i) {
1358                            o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;
1359                            o[INDEX2(i,numComp,1)] = f_01[i]*c0 + f_11[i]*c1;
1360                        } /* end of component loop i */
1361                    } /* end of k0 loop */
1362                } /* end of face 3 */
1363                /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
1364            } // end of parallel section
1365      }      }
1366    }
1367    
1368      cout << "--- colCouple_pattern ---" << endl;  //protected
1369      cout << "M=" << M << ", N=" << N << endl;  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1370      for (size_t i=0; i<ptr.size(); i++) {          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1371          cout << "ptr[" << i << "]=" << ptr[i] << endl;          const escript::Data& C, const escript::Data& D,
1372      }          const escript::Data& X, const escript::Data& Y,
1373      for (size_t i=0; i<index.size(); i++) {          const escript::Data& d, const escript::Data& y) const
1374          cout << "index[" << i << "]=" << index[i] << endl;  {
1375      }      const double h0 = m_l0/m_gNE0;
1376        const double h1 = m_l1/m_gNE1;
1377      // now build the row couple pattern      /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */
1378      IndexVector ptr2(1,0);      const double w0 = -0.1555021169820365539*h1/h0;
1379      IndexVector index2;      const double w1 = 0.041666666666666666667;
1380      for (dim_t id=0; id<N; id++) {      const double w10 = -0.041666666666666666667*h0/h1;
1381          n=0;      const double w11 = 0.1555021169820365539*h1/h0;
1382          for (dim_t i=0; i<M; i++) {      const double w12 = 0.1555021169820365539*h0/h1;
1383              for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {      const double w13 = 0.01116454968463011277*h0/h1;
1384                  if (index[j]==id) {      const double w14 = 0.01116454968463011277*h1/h0;
1385                      index2.push_back(i);      const double w15 = 0.041666666666666666667*h1/h0;
1386                      n++;      const double w16 = -0.01116454968463011277*h0/h1;
1387                      break;      const double w17 = -0.1555021169820365539*h0/h1;
1388                  }      const double w18 = -0.33333333333333333333*h1/h0;
1389              }      const double w19 = 0.25000000000000000000;
1390          }      const double w2 = -0.15550211698203655390;
1391          ptr2.push_back(ptr2.back()+n);      const double w20 = -0.25000000000000000000;
1392      }      const double w21 = 0.16666666666666666667*h0/h1;
1393        const double w22 = -0.16666666666666666667*h1/h0;
1394        const double w23 = -0.16666666666666666667*h0/h1;
1395        const double w24 = 0.33333333333333333333*h1/h0;
1396        const double w25 = 0.33333333333333333333*h0/h1;
1397        const double w26 = 0.16666666666666666667*h1/h0;
1398        const double w27 = -0.33333333333333333333*h0/h1;
1399        const double w28 = -0.032861463941450536761*h1;
1400        const double w29 = -0.032861463941450536761*h0;
1401        const double w3 = 0.041666666666666666667*h0/h1;
1402        const double w30 = -0.12264065304058601714*h1;
1403        const double w31 = -0.0023593469594139828636*h1;
1404        const double w32 = -0.008805202725216129906*h0;
1405        const double w33 = -0.008805202725216129906*h1;
1406        const double w34 = 0.032861463941450536761*h1;
1407        const double w35 = 0.008805202725216129906*h1;
1408        const double w36 = 0.008805202725216129906*h0;
1409        const double w37 = 0.0023593469594139828636*h1;
1410        const double w38 = 0.12264065304058601714*h1;
1411        const double w39 = 0.032861463941450536761*h0;
1412        const double w4 = 0.15550211698203655390;
1413        const double w40 = -0.12264065304058601714*h0;
1414        const double w41 = -0.0023593469594139828636*h0;
1415        const double w42 = 0.0023593469594139828636*h0;
1416        const double w43 = 0.12264065304058601714*h0;
1417        const double w44 = -0.16666666666666666667*h1;
1418        const double w45 = -0.083333333333333333333*h0;
1419        const double w46 = 0.083333333333333333333*h1;
1420        const double w47 = 0.16666666666666666667*h1;
1421        const double w48 = 0.083333333333333333333*h0;
1422        const double w49 = -0.16666666666666666667*h0;
1423        const double w5 = -0.041666666666666666667;
1424        const double w50 = 0.16666666666666666667*h0;
1425        const double w51 = -0.083333333333333333333*h1;
1426        const double w52 = 0.025917019497006092316*h0*h1;
1427        const double w53 = 0.0018607582807716854616*h0*h1;
1428        const double w54 = 0.0069444444444444444444*h0*h1;
1429        const double w55 = 0.09672363354357992482*h0*h1;
1430        const double w56 = 0.00049858867864229740201*h0*h1;
1431        const double w57 = 0.055555555555555555556*h0*h1;
1432        const double w58 = 0.027777777777777777778*h0*h1;
1433        const double w59 = 0.11111111111111111111*h0*h1;
1434        const double w6 = -0.01116454968463011277*h1/h0;
1435        const double w60 = -0.19716878364870322056*h1;
1436        const double w61 = -0.19716878364870322056*h0;
1437        const double w62 = -0.052831216351296779436*h0;
1438        const double w63 = -0.052831216351296779436*h1;
1439        const double w64 = 0.19716878364870322056*h1;
1440        const double w65 = 0.052831216351296779436*h1;
1441        const double w66 = 0.19716878364870322056*h0;
1442        const double w67 = 0.052831216351296779436*h0;
1443        const double w68 = -0.5*h1;
1444        const double w69 = -0.5*h0;
1445        const double w7 = 0.011164549684630112770;
1446        const double w70 = 0.5*h1;
1447        const double w71 = 0.5*h0;
1448        const double w72 = 0.1555021169820365539*h0*h1;
1449        const double w73 = 0.041666666666666666667*h0*h1;
1450        const double w74 = 0.01116454968463011277*h0*h1;
1451        const double w75 = 0.25*h0*h1;
1452        const double w8 = -0.011164549684630112770;
1453        const double w9 = -0.041666666666666666667*h1/h0;
1454        /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */
1455    
1456      cout << "--- rowCouple_pattern ---" << endl;      rhs.requireWrite();
1457      cout << "M=" << N << ", N=" << M << endl;  #pragma omp parallel
1458      for (size_t i=0; i<ptr2.size(); i++) {      {
1459          cout << "ptr[" << i << "]=" << ptr2[i] << endl;          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1460      }  #pragma omp for
1461      for (size_t i=0; i<index2.size(); i++) {              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
1462          cout << "index[" << i << "]=" << index2[i] << endl;                  for (index_t k0=0; k0<m_NE0; ++k0)  {
1463      }                      bool add_EM_S=false;
1464                        bool add_EM_F=false;
1465      // paso will manage the memory                      vector<double> EM_S(4*4, 0);
1466      index_t* indexC = MEMALLOC(index.size(), index_t);                      vector<double> EM_F(4, 0);
1467      index_t* ptrC = MEMALLOC(ptr.size(), index_t);                      const index_t e = k0 + m_NE0*k1;
1468      copy(index.begin(), index.end(), indexC);                      /*** GENERATOR SNIP_PDE_SINGLE TOP */
1469      copy(ptr.begin(), ptr.end(), ptrC);                      ///////////////
1470      *colPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);                      // process A //
1471                        ///////////////
1472      // paso will manage the memory                      if (!A.isEmpty()) {
1473      indexC = MEMALLOC(index2.size(), index_t);                          add_EM_S=true;
1474      ptrC = MEMALLOC(ptr2.size(), index_t);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1475      copy(index2.begin(), index2.end(), indexC);                          if (A.actsExpanded()) {
1476      copy(ptr2.begin(), ptr2.end(), ptrC);                              const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1477      *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);                              const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1478                                const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1479                                const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1480                                const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1481                                const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1482                                const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1483                                const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1484                                const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1485                                const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1486                                const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1487                                const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1488                                const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1489                                const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1490                                const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1491                                const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1492                                const register double tmp4_0 = A_10_1 + A_10_2;
1493                                const register double tmp12_0 = A_11_0 + A_11_2;
1494                                const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1495                                const register double tmp10_0 = A_01_3 + A_10_3;
1496                                const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
1497                                const register double tmp0_0 = A_01_0 + A_01_3;
1498                                const register double tmp13_0 = A_01_2 + A_10_1;
1499                                const register double tmp3_0 = A_00_2 + A_00_3;
1500                                const register double tmp11_0 = A_11_1 + A_11_3;
1501                                const register double tmp18_0 = A_01_1 + A_10_1;
1502                                const register double tmp1_0 = A_00_0 + A_00_1;
1503                                const register double tmp15_0 = A_01_1 + A_10_2;
1504                                const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1505                                const register double tmp16_0 = A_10_0 + A_10_3;
1506                                const register double tmp6_0 = A_01_3 + A_10_0;
1507                                const register double tmp17_0 = A_01_1 + A_01_2;
1508                                const register double tmp9_0 = A_01_0 + A_10_0;
1509                                const register double tmp7_0 = A_01_0 + A_10_3;
1510                                const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1511                                const register double tmp19_0 = A_01_2 + A_10_2;
1512                                const register double tmp14_1 = A_10_0*w8;
1513                                const register double tmp23_1 = tmp3_0*w14;
1514                                const register double tmp35_1 = A_01_0*w8;
1515                                const register double tmp54_1 = tmp13_0*w8;
1516                                const register double tmp20_1 = tmp9_0*w4;
1517                                const register double tmp25_1 = tmp12_0*w12;
1518                                const register double tmp2_1 = A_01_1*w4;
1519                                const register double tmp44_1 = tmp7_0*w7;
1520                                const register double tmp26_1 = tmp10_0*w4;
1521                                const register double tmp52_1 = tmp18_0*w8;
1522                                const register double tmp48_1 = A_10_1*w7;
1523                                const register double tmp46_1 = A_01_3*w8;
1524                                const register double tmp50_1 = A_01_0*w2;
1525                                const register double tmp8_1 = tmp4_0*w5;
1526                                const register double tmp56_1 = tmp19_0*w8;
1527                                const register double tmp9_1 = tmp2_0*w10;
1528                                const register double tmp19_1 = A_10_3*w2;
1529                                const register double tmp47_1 = A_10_2*w4;
1530                                const register double tmp16_1 = tmp3_0*w0;
1531                                const register double tmp18_1 = tmp1_0*w6;
1532                                const register double tmp31_1 = tmp11_0*w12;
1533                                const register double tmp55_1 = tmp15_0*w2;
1534                                const register double tmp39_1 = A_10_2*w7;
1535                                const register double tmp11_1 = tmp6_0*w7;
1536                                const register double tmp40_1 = tmp11_0*w17;
1537                                const register double tmp34_1 = tmp15_0*w8;
1538                                const register double tmp33_1 = tmp14_0*w5;
1539                                const register double tmp24_1 = tmp11_0*w13;
1540                                const register double tmp3_1 = tmp1_0*w0;
1541                                const register double tmp5_1 = tmp2_0*w3;
1542                                const register double tmp43_1 = tmp17_0*w5;
1543                                const register double tmp15_1 = A_01_2*w4;
1544                                const register double tmp53_1 = tmp19_0*w2;
1545                                const register double tmp27_1 = tmp3_0*w11;
1546                                const register double tmp32_1 = tmp13_0*w2;
1547                                const register double tmp10_1 = tmp5_0*w9;
1548                                const register double tmp37_1 = A_10_1*w4;
1549                                const register double tmp38_1 = tmp5_0*w15;
1550                                const register double tmp17_1 = A_01_1*w7;
1551                                const register double tmp12_1 = tmp7_0*w4;
1552                                const register double tmp22_1 = tmp10_0*w7;
1553                                const register double tmp57_1 = tmp18_0*w2;
1554                                const register double tmp28_1 = tmp9_0*w7;
1555                                const register double tmp29_1 = tmp1_0*w14;
1556                                const register double tmp51_1 = tmp11_0*w16;
1557                                const register double tmp42_1 = tmp12_0*w16;
1558                                const register double tmp49_1 = tmp12_0*w17;
1559                                const register double tmp21_1 = tmp1_0*w11;
1560                                const register double tmp1_1 = tmp0_0*w1;
1561                                const register double tmp45_1 = tmp6_0*w4;
1562                                const register double tmp7_1 = A_10_0*w2;
1563                                const register double tmp6_1 = tmp3_0*w6;
1564                                const register double tmp13_1 = tmp8_0*w1;
1565                                const register double tmp36_1 = tmp16_0*w1;
1566                                const register double tmp41_1 = A_01_3*w2;
1567                                const register double tmp30_1 = tmp12_0*w13;
1568                                const register double tmp4_1 = A_01_2*w7;
1569                                const register double tmp0_1 = A_10_3*w8;
1570                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
1571                                EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
1572                                EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
1573                                EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1574                                EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1575                                EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
1576                                EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
1577                                EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
1578                                EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1579                                EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
1580                                EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
1581                                EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
1582                                EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
1583                                EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
1584                                EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
1585                                EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
1586                            } else { /* constant data */
1587                                const register double A_00 = A_p[INDEX2(0,0,2)];
1588                                const register double A_01 = A_p[INDEX2(0,1,2)];
1589                                const register double A_10 = A_p[INDEX2(1,0,2)];
1590                                const register double A_11 = A_p[INDEX2(1,1,2)];
1591                                const register double tmp0_0 = A_01 + A_10;
1592                                const register double tmp0_1 = A_00*w18;
1593                                const register double tmp10_1 = A_01*w20;
1594                                const register double tmp12_1 = A_00*w26;
1595                                const register double tmp4_1 = A_00*w22;
1596                                const register double tmp8_1 = A_00*w24;
1597                                const register double tmp13_1 = A_10*w19;
1598                                const register double tmp9_1 = tmp0_0*w20;
1599                                const register double tmp3_1 = A_11*w21;
1600                                const register double tmp11_1 = A_11*w27;
1601                                const register double tmp1_1 = A_01*w19;
1602                                const register double tmp6_1 = A_11*w23;
1603                                const register double tmp7_1 = A_11*w25;
1604                                const register double tmp2_1 = A_10*w20;
1605                                const register double tmp5_1 = tmp0_0*w19;
1606                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1607                                EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1608                                EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1609                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1610                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1611                                EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1612                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1613                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1614                                EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1615                                EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1616                                EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1617                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1618                                EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1619                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1620                                EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1621                                EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1622                            }
1623                        }
1624                        ///////////////
1625                        // process B //
1626                        ///////////////
1627                        if (!B.isEmpty()) {
1628                            add_EM_S=true;
1629                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1630                            if (B.actsExpanded()) {
1631                                const register double B_0_0 = B_p[INDEX2(0,0,2)];
1632                                const register double B_1_0 = B_p[INDEX2(1,0,2)];
1633                                const register double B_0_1 = B_p[INDEX2(0,1,2)];
1634                                const register double B_1_1 = B_p[INDEX2(1,1,2)];
1635                                const register double B_0_2 = B_p[INDEX2(0,2,2)];
1636                                const register double B_1_2 = B_p[INDEX2(1,2,2)];
1637                                const register double B_0_3 = B_p[INDEX2(0,3,2)];
1638                                const register double B_1_3 = B_p[INDEX2(1,3,2)];
1639                                const register double tmp3_0 = B_0_0 + B_0_2;
1640                                const register double tmp1_0 = B_1_2 + B_1_3;
1641                                const register double tmp2_0 = B_0_1 + B_0_3;
1642                                const register double tmp0_0 = B_1_0 + B_1_1;
1643                                const register double tmp63_1 = B_1_1*w42;
1644                                const register double tmp79_1 = B_1_1*w40;
1645                                const register double tmp37_1 = tmp3_0*w35;
1646                                const register double tmp8_1 = tmp0_0*w32;
1647                                const register double tmp71_1 = B_0_1*w34;
1648                                const register double tmp19_1 = B_0_3*w31;
1649                                const register double tmp15_1 = B_0_3*w34;
1650                                const register double tmp9_1 = tmp3_0*w34;
1651                                const register double tmp35_1 = B_1_0*w36;
1652                                const register double tmp66_1 = B_0_3*w28;
1653                                const register double tmp28_1 = B_1_0*w42;
1654                                const register double tmp22_1 = B_1_0*w40;
1655                                const register double tmp16_1 = B_1_2*w29;
1656                                const register double tmp6_1 = tmp2_0*w35;
1657                                const register double tmp55_1 = B_1_3*w40;
1658                                const register double tmp50_1 = B_1_3*w42;
1659                                const register double tmp7_1 = tmp1_0*w29;
1660                                const register double tmp1_1 = tmp1_0*w32;
1661                                const register double tmp57_1 = B_0_3*w30;
1662                                const register double tmp18_1 = B_1_1*w32;
1663                                const register double tmp53_1 = B_1_0*w41;
1664                                const register double tmp61_1 = B_1_3*w36;
1665                                const register double tmp27_1 = B_0_3*w38;
1666                                const register double tmp64_1 = B_0_2*w30;
1667                                const register double tmp76_1 = B_0_1*w38;
1668                                const register double tmp39_1 = tmp2_0*w34;
1669                                const register double tmp62_1 = B_0_1*w31;
1670                                const register double tmp56_1 = B_0_0*w31;
1671                                const register double tmp49_1 = B_1_1*w36;
1672                                const register double tmp2_1 = B_0_2*w31;
1673                                const register double tmp23_1 = B_0_2*w33;
1674                                const register double tmp38_1 = B_1_1*w43;
1675                                const register double tmp74_1 = B_1_2*w41;
1676                                const register double tmp43_1 = B_1_1*w41;
1677                                const register double tmp58_1 = B_0_2*w28;
1678                                const register double tmp67_1 = B_0_0*w33;
1679                                const register double tmp33_1 = tmp0_0*w39;
1680                                const register double tmp4_1 = B_0_0*w28;
1681                                const register double tmp20_1 = B_0_0*w30;
1682                                const register double tmp13_1 = B_0_2*w38;
1683                                const register double tmp65_1 = B_1_2*w43;
1684                                const register double tmp0_1 = tmp0_0*w29;
1685                                const register double tmp41_1 = tmp3_0*w33;
1686                                const register double tmp73_1 = B_0_2*w37;
1687                                const register double tmp69_1 = B_0_0*w38;
1688                                const register double tmp48_1 = B_1_2*w39;
1689                                const register double tmp59_1 = B_0_1*w33;
1690                                const register double tmp17_1 = B_1_3*w41;
1691                                const register double tmp5_1 = B_0_3*w33;
1692                                const register double tmp3_1 = B_0_1*w30;
1693                                const register double tmp21_1 = B_0_1*w28;
1694                                const register double tmp42_1 = B_1_0*w29;
1695                                const register double tmp54_1 = B_1_2*w32;
1696                                const register double tmp60_1 = B_1_0*w39;
1697                                const register double tmp32_1 = tmp1_0*w36;
1698                                const register double tmp10_1 = B_0_1*w37;
1699                                const register double tmp14_1 = B_0_0*w35;
1700                                const register double tmp29_1 = B_0_1*w35;
1701                                const register double tmp26_1 = B_1_2*w36;
1702                                const register double tmp30_1 = B_1_3*w43;
1703                                const register double tmp70_1 = B_0_2*w35;
1704                                const register double tmp34_1 = B_1_3*w39;
1705                                const register double tmp51_1 = B_1_0*w43;
1706                                const register double tmp31_1 = B_0_2*w34;
1707                                const register double tmp45_1 = tmp3_0*w28;
1708                                const register double tmp11_1 = tmp1_0*w39;
1709                                const register double tmp52_1 = B_1_1*w29;
1710                                const register double tmp44_1 = B_1_3*w32;
1711                                const register double tmp25_1 = B_1_1*w39;
1712                                const register double tmp47_1 = tmp2_0*w33;
1713                                const register double tmp72_1 = B_1_3*w29;
1714                                const register double tmp40_1 = tmp2_0*w28;
1715                                const register double tmp46_1 = B_1_2*w40;
1716                                const register double tmp36_1 = B_1_2*w42;
1717                                const register double tmp24_1 = B_0_0*w37;
1718                                const register double tmp77_1 = B_0_3*w35;
1719                                const register double tmp68_1 = B_0_3*w37;
1720                                const register double tmp78_1 = B_0_0*w34;
1721                                const register double tmp12_1 = tmp0_0*w36;
1722                                const register double tmp75_1 = B_1_0*w32;
1723                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1724                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1725                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1726                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1727                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1728                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
1729                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1730                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1731                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1732                                EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1733                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1734                                EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1735                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1736                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1737                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
1738                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1739                            } else { /* constant data */
1740                                const register double B_0 = B_p[0];
1741                                const register double B_1 = B_p[1];
1742                                const register double tmp6_1 = B_1*w50;
1743                                const register double tmp1_1 = B_1*w45;
1744                                const register double tmp5_1 = B_1*w49;
1745                                const register double tmp4_1 = B_1*w48;
1746                                const register double tmp0_1 = B_0*w44;
1747                                const register double tmp2_1 = B_0*w46;
1748                                const register double tmp7_1 = B_0*w51;
1749                                const register double tmp3_1 = B_0*w47;
1750                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1751                                EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
1752                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1753                                EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1754                                EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
1755                                EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
1756                                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
1757                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
1758                                EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
1759                                EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
1760                                EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
1761                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
1762                                EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
1763                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1764                                EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
1765                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
1766                            }
1767                        }
1768                        ///////////////
1769                        // process C //
1770                        ///////////////
1771                        if (!C.isEmpty()) {
1772                            add_EM_S=true;
1773                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
1774                            if (C.actsExpanded()) {
1775                                const register double C_0_0 = C_p[INDEX2(0,0,2)];
1776                                const register double C_1_0 = C_p[INDEX2(1,0,2)];
1777                                const register double C_0_1 = C_p[INDEX2(0,1,2)];
1778                                const register double C_1_1 = C_p[INDEX2(1,1,2)];
1779                                const register double C_0_2 = C_p[INDEX2(0,2,2)];
1780                                const register double C_1_2 = C_p[INDEX2(1,2,2)];
1781                                const register double C_0_3 = C_p[INDEX2(0,3,2)];
1782                                const register double C_1_3 = C_p[INDEX2(1,3,2)];
1783                                const register double tmp2_0 = C_0_1 + C_0_3;
1784                                const register double tmp1_0 = C_1_2 + C_1_3;
1785                                const register double tmp3_0 = C_0_0 + C_0_2;
1786                                const register double tmp0_0 = C_1_0 + C_1_1;
1787                                const register double tmp64_1 = C_0_2*w30;
1788                                const register double tmp14_1 = C_0_2*w28;
1789                                const register double tmp19_1 = C_0_3*w31;
1790                                const register double tmp22_1 = C_1_0*w40;
1791                                const register double tmp37_1 = tmp3_0*w35;
1792                                const register double tmp29_1 = C_0_1*w35;
1793                                const register double tmp73_1 = C_0_2*w37;
1794                                const register double tmp74_1 = C_1_2*w41;
1795                                const register double tmp52_1 = C_1_3*w39;
1796                                const register double tmp25_1 = C_1_1*w39;
1797                                const register double tmp62_1 = C_0_1*w31;
1798                                const register double tmp79_1 = C_1_1*w40;
1799                                const register double tmp43_1 = C_1_1*w36;
1800                                const register double tmp27_1 = C_0_3*w38;
1801                                const register double tmp28_1 = C_1_0*w42;
1802                                const register double tmp63_1 = C_1_1*w42;
1803                                const register double tmp59_1 = C_0_3*w34;
1804                                const register double tmp72_1 = C_1_3*w29;
1805                                const register double tmp40_1 = tmp2_0*w35;
1806                                const register double tmp13_1 = C_0_3*w30;
1807                                const register double tmp51_1 = C_1_2*w40;
1808                                const register double tmp54_1 = C_1_2*w42;
1809                                const register double tmp12_1 = C_0_0*w31;
1810                                const register double tmp2_1 = tmp1_0*w32;
1811                                const register double tmp68_1 = C_0_2*w31;
1812                                const register double tmp75_1 = C_1_0*w32;
1813                                const register double tmp49_1 = C_1_1*w41;
1814                                const register double tmp4_1 = C_0_2*w35;
1815                                const register double tmp66_1 = C_0_3*w28;
1816                                const register double tmp56_1 = C_0_1*w37;
1817                                const register double tmp5_1 = C_0_1*w34;
1818                                const register double tmp38_1 = tmp2_0*w34;
1819                                const register double tmp76_1 = C_0_1*w38;
1820                                const register double tmp21_1 = C_0_1*w28;
1821                                const register double tmp69_1 = C_0_1*w30;
1822                                const register double tmp53_1 = C_1_0*w36;
1823                                const register double tmp42_1 = C_1_2*w39;
1824                                const register double tmp32_1 = tmp1_0*w29;
1825                                const register double tmp45_1 = C_1_0*w43;
1826                                const register double tmp33_1 = tmp0_0*w32;
1827                                const register double tmp35_1 = C_1_0*w41;
1828                                const register double tmp26_1 = C_1_2*w36;
1829                                const register double tmp67_1 = C_0_0*w33;
1830                                const register double tmp31_1 = C_0_2*w34;
1831                                const register double tmp20_1 = C_0_0*w30;
1832                                const register double tmp70_1 = C_0_0*w28;
1833                                const register double tmp8_1 = tmp0_0*w39;
1834                                const register double tmp30_1 = C_1_3*w43;
1835                                const register double tmp0_1 = tmp0_0*w29;
1836                                const register double tmp17_1 = C_1_3*w41;
1837                                const register double tmp58_1 = C_0_0*w35;
1838                                const register double tmp9_1 = tmp3_0*w33;
1839                                const register double tmp61_1 = C_1_3*w36;
1840                                const register double tmp41_1 = tmp3_0*w34;
1841                                const register double tmp50_1 = C_1_3*w32;
1842                                const register double tmp18_1 = C_1_1*w32;
1843                                const register double tmp6_1 = tmp1_0*w36;
1844                                const register double tmp3_1 = C_0_0*w38;
1845                                const register double tmp34_1 = C_1_1*w29;
1846                                const register double tmp77_1 = C_0_3*w35;
1847                                const register double tmp65_1 = C_1_2*w43;
1848                                const register double tmp71_1 = C_0_3*w33;
1849                                const register double tmp55_1 = C_1_1*w43;
1850                                const register double tmp46_1 = tmp3_0*w28;
1851                                const register double tmp24_1 = C_0_0*w37;
1852                                const register double tmp10_1 = tmp1_0*w39;
1853                                const register double tmp48_1 = C_1_0*w29;
1854                                const register double tmp15_1 = C_0_1*w33;
1855                                const register double tmp36_1 = C_1_2*w32;
1856                                const register double tmp60_1 = C_1_0*w39;
1857                                const register double tmp47_1 = tmp2_0*w33;
1858                                const register double tmp16_1 = C_1_2*w29;
1859                                const register double tmp1_1 = C_0_3*w37;
1860                                const register double tmp7_1 = tmp2_0*w28;
1861                                const register double tmp39_1 = C_1_3*w40;
1862                                const register double tmp44_1 = C_1_3*w42;
1863                                const register double tmp57_1 = C_0_2*w38;
1864                                const register double tmp78_1 = C_0_0*w34;
1865                                const register double tmp11_1 = tmp0_0*w36;
1866                                const register double tmp23_1 = C_0_2*w33;
1867                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1868                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1869                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1870                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1871                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1872                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
1873                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1874                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1875                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1876                                EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1877                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1878                                EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1879                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1880                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1881                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
1882                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1883                            } else { /* constant data */
1884                                const register double C_0 = C_p[0];
1885                                const register double C_1 = C_p[1];
1886                                const register double tmp1_1 = C_1*w45;
1887                                const register double tmp3_1 = C_0*w51;
1888                                const register double tmp4_1 = C_0*w44;
1889                                const register double tmp7_1 = C_0*w46;
1890                                const register double tmp5_1 = C_1*w49;
1891                                const register double tmp2_1 = C_1*w48;
1892                                const register double tmp0_1 = C_0*w47;
1893                                const register double tmp6_1 = C_1*w50;
1894                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1895                                EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
1896                                EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
1897                                EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1898                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
1899                                EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
1900                                EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
1901                                EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
1902                                EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
1903                                EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
1904                                EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
1905                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
1906                                EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
1907                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
1908                                EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
1909                                EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
1910                            }
1911                        }
1912                        ///////////////
1913                        // process D //
1914                        ///////////////
1915                        if (!D.isEmpty()) {
1916                            add_EM_S=true;
1917                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
1918                            if (D.actsExpanded()) {
1919                                const register double D_0 = D_p[0];
1920                                const register double D_1 = D_p[1];
1921                                const register double D_2 = D_p[2];
1922                                const register double D_3 = D_p[3];
1923                                const register double tmp4_0 = D_1 + D_3;
1924                                const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;
1925                                const register double tmp5_0 = D_0 + D_2;
1926                                const register double tmp0_0 = D_0 + D_1;
1927                                const register double tmp6_0 = D_0 + D_3;
1928                                const register double tmp1_0 = D_2 + D_3;
1929                                const register double tmp3_0 = D_1 + D_2;
1930                                const register double tmp16_1 = D_1*w56;
1931                                const register double tmp14_1 = tmp6_0*w54;
1932                                const register double tmp8_1 = D_3*w55;
1933                                const register double tmp2_1 = tmp2_0*w54;
1934                                const register double tmp12_1 = tmp5_0*w52;
1935                                const register double tmp4_1 = tmp0_0*w53;
1936                                const register double tmp3_1 = tmp1_0*w52;
1937                                const register double tmp13_1 = tmp4_0*w53;
1938                                const register double tmp10_1 = tmp4_0*w52;
1939                                const register double tmp1_1 = tmp1_0*w53;
1940                                const register double tmp7_1 = D_3*w56;
1941                                const register double tmp0_1 = tmp0_0*w52;
1942                                const register double tmp11_1 = tmp5_0*w53;
1943                                const register double tmp9_1 = D_0*w56;
1944                                const register double tmp5_1 = tmp3_0*w54;
1945                                const register double tmp18_1 = D_2*w56;
1946                                const register double tmp17_1 = D_1*w55;
1947                                const register double tmp6_1 = D_0*w55;
1948                                const register double tmp15_1 = D_2*w55;
1949                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1950                                EM_S[INDEX2(1,2,4)]+=tmp2_1;
1951                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1952                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
1953                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
1954                                EM_S[INDEX2(3,0,4)]+=tmp2_1;
1955                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
1956                                EM_S[INDEX2(2,1,4)]+=tmp2_1;
1957                                EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
1958                                EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
1959                                EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
1960                                EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
1961                                EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
1962                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
1963                                EM_S[INDEX2(0,3,4)]+=tmp2_1;
1964                                EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
1965                            } else { /* constant data */
1966                                const register double D_0 = D_p[0];
1967                                const register double tmp2_1 = D_0*w59;
1968                                const register double tmp1_1 = D_0*w58;
1969                                const register double tmp0_1 = D_0*w57;
1970                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
1971                                EM_S[INDEX2(1,2,4)]+=tmp1_1;
1972                                EM_S[INDEX2(3,2,4)]+=tmp0_1;
1973                                EM_S[INDEX2(0,0,4)]+=tmp2_1;
1974                                EM_S[INDEX2(3,3,4)]+=tmp2_1;
1975                                EM_S[INDEX2(3,0,4)]+=tmp1_1;
1976                                EM_S[INDEX2(3,1,4)]+=tmp0_1;
1977                                EM_S[INDEX2(2,1,4)]+=tmp1_1;
1978                                EM_S[INDEX2(0,2,4)]+=tmp0_1;
1979                                EM_S[INDEX2(2,0,4)]+=tmp0_1;
1980                                EM_S[INDEX2(1,3,4)]+=tmp0_1;
1981                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
1982                                EM_S[INDEX2(2,2,4)]+=tmp2_1;
1983                                EM_S[INDEX2(1,0,4)]+=tmp0_1;
1984                                EM_S[INDEX2(0,3,4)]+=tmp1_1;
1985                                EM_S[INDEX2(1,1,4)]+=tmp2_1;
1986                            }
1987                        }
1988                        ///////////////
1989                        // process X //
1990                        ///////////////
1991                        if (!X.isEmpty()) {
1992                            add_EM_F=true;
1993                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
1994                            if (X.actsExpanded()) {
1995                                const register double X_0_0 = X_p[INDEX2(0,0,2)];
1996                                const register double X_1_0 = X_p[INDEX2(1,0,2)];
1997                                const register double X_0_1 = X_p[INDEX2(0,1,2)];
1998                                const register double X_1_1 = X_p[INDEX2(1,1,2)];
1999                                const register double X_0_2 = X_p[INDEX2(0,2,2)];
2000                                const register double X_1_2 = X_p[INDEX2(1,2,2)];
2001                                const register double X_0_3 = X_p[INDEX2(0,3,2)];
2002                                const register double X_1_3 = X_p[INDEX2(1,3,2)];
2003                                const register double tmp1_0 = X_1_1 + X_1_3;
2004                                const register double tmp3_0 = X_0_0 + X_0_1;
2005                                const register double tmp2_0 = X_1_0 + X_1_2;
2006                                const register double tmp0_0 = X_0_2 + X_0_3;
2007                                const register double tmp8_1 = tmp2_0*w66;
2008                                const register double tmp5_1 = tmp3_0*w64;
2009                                const register double tmp14_1 = tmp0_0*w64;
2010                                const register double tmp3_1 = tmp3_0*w60;
2011                                const register double tmp9_1 = tmp3_0*w63;
2012                                const register double tmp13_1 = tmp3_0*w65;
2013                                const register double tmp12_1 = tmp1_0*w66;
2014                                const register double tmp10_1 = tmp0_0*w60;
2015                                const register double tmp2_1 = tmp2_0*w61;
2016                                const register double tmp6_1 = tmp2_0*w62;
2017                                const register double tmp4_1 = tmp0_0*w65;
2018                                const register double tmp11_1 = tmp1_0*w67;
2019                                const register double tmp1_1 = tmp1_0*w62;
2020                                const register double tmp7_1 = tmp1_0*w61;
2021                                const register double tmp0_1 = tmp0_0*w63;
2022                                const register double tmp15_1 = tmp2_0*w67;
2023                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2024                                EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2025                                EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2026                                EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2027                            } else { /* constant data */
2028                                const register double X_0 = X_p[0];
2029                                const register double X_1 = X_p[1];
2030                                const register double tmp3_1 = X_1*w71;
2031                                const register double tmp2_1 = X_0*w70;
2032                                const register double tmp1_1 = X_0*w68;
2033                                const register double tmp0_1 = X_1*w69;
2034                                EM_F[0]+=tmp0_1 + tmp1_1;
2035                                EM_F[1]+=tmp0_1 + tmp2_1;
2036                                EM_F[2]+=tmp1_1 + tmp3_1;
2037                                EM_F[3]+=tmp2_1 + tmp3_1;
2038                            }
2039                        }
2040                        ///////////////
2041                        // process Y //
2042                        ///////////////
2043                        if (!Y.isEmpty()) {
2044                            add_EM_F=true;
2045                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2046                            if (Y.actsExpanded()) {
2047                                const register double Y_0 = Y_p[0];
2048                                const register double Y_1 = Y_p[1];
2049                                const register double Y_2 = Y_p[2];
2050                                const register double Y_3 = Y_p[3];
2051                                const register double tmp0_0 = Y_1 + Y_2;
2052                                const register double tmp1_0 = Y_0 + Y_3;
2053                                const register double tmp9_1 = Y_0*w74;
2054                                const register double tmp4_1 = tmp1_0*w73;
2055                                const register double tmp5_1 = Y_2*w74;
2056                                const register double tmp7_1 = Y_1*w74;
2057                                const register double tmp6_1 = Y_2*w72;
2058                                const register double tmp2_1 = Y_3*w74;
2059                                const register double tmp8_1 = Y_3*w72;
2060                                const register double tmp3_1 = Y_1*w72;
2061                                const register double tmp0_1 = Y_0*w72;
2062                                const register double tmp1_1 = tmp0_0*w73;
2063                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2064                                EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2065                                EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2066                                EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2067                            } else { /* constant data */
2068                                const register double Y_0 = Y_p[0];
2069                                const register double tmp0_1 = Y_0*w75;
2070                                EM_F[0]+=tmp0_1;
2071                                EM_F[1]+=tmp0_1;
2072                                EM_F[2]+=tmp0_1;
2073                                EM_F[3]+=tmp0_1;
2074                            }
2075                        }
2076                        /* GENERATOR SNIP_PDE_SINGLE BOTTOM */
2077    
2078                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2079                        const index_t firstNode=m_N0*k1+k0;
2080                        IndexVector rowIndex;
2081                        rowIndex.push_back(m_dofMap[firstNode]);
2082                        rowIndex.push_back(m_dofMap[firstNode+1]);
2083                        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2084                        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2085                        if (add_EM_F) {
2086                            //cout << "-- AddtoRHS -- " << endl;
2087                            double *F_p=rhs.getSampleDataRW(0);
2088                            for (index_t i=0; i<rowIndex.size(); i++) {
2089                                if (rowIndex[i]<getNumDOF()) {
2090                                    F_p[rowIndex[i]]+=EM_F[i];
2091                                    //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;
2092                                }
2093                            }
2094                            //cout << "---"<<endl;
2095                        }
2096                        if (add_EM_S) {
2097                            //cout << "-- AddtoSystem -- " << endl;
2098                            addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2099                        }
2100    
2101                    } // end k0 loop
2102                } // end k1 loop
2103            } // end of coloring
2104        } // end of parallel region
2105  }  }
2106    
2107  } // end of namespace ripley  } // end of namespace ripley

Legend:
Removed from v.3699  
changed lines
  Added in v.3759

  ViewVC Help
Powered by ViewVC 1.1.26