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

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

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

revision 3697 by caltinay, Mon Nov 28 04:52:00 2011 UTC revision 3748 by caltinay, Thu Dec 15 07:36:19 2011 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 69  string Rectangle::getDescription() const Line 69  string Rectangle::getDescription() const
69    
70  bool Rectangle::operator==(const AbstractDomain& other) const  bool Rectangle::operator==(const AbstractDomain& other) const
71  {  {
72      if (dynamic_cast<const Rectangle*>(&other))      const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
73          return this==&other;      if (o) {
74            return (RipleyDomain::operator==(other) &&
75                    m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1
76                    && m_l0==o->m_l0 && m_l1==o->m_l1
77                    && m_NX==o->m_NX && m_NY==o->m_NY);
78        }
79    
80      return false;      return false;
81  }  }
# Line 142  void Rectangle::dump(const string& fileN Line 147  void Rectangle::dump(const string& fileN
147      pair<double,double> ydy = getFirstCoordAndSpacing(1);      pair<double,double> ydy = getFirstCoordAndSpacing(1);
148  #pragma omp parallel  #pragma omp parallel
149      {      {
150  #pragma omp for  #pragma omp for nowait
151          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_N0; i0++) {
152              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=xdx.first+i0*xdx.second;
153          }          }
154  #pragma omp for  #pragma omp for nowait
155          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_N1; i1++) {
156              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=ydy.first+i1*ydy.second;
157          }          }
# Line 221  const int* Rectangle::borrowSampleRefere Line 226  const int* Rectangle::borrowSampleRefere
226  {  {
227      switch (fsType) {      switch (fsType) {
228          case Nodes:          case Nodes:
229            case ReducedNodes: //FIXME: reduced
230              return &m_nodeId[0];              return &m_nodeId[0];
231          case Elements:          case Elements:
232            case ReducedElements:
233              return &m_elementId[0];              return &m_elementId[0];
234          case FaceElements:          case FaceElements:
235            case ReducedFaceElements:
236              return &m_faceId[0];              return &m_faceId[0];
237          default:          default:
238              break;              break;
# Line 232  const int* Rectangle::borrowSampleRefere Line 240  const int* Rectangle::borrowSampleRefere
240    
241      stringstream msg;      stringstream msg;
242      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs() not implemented for function space type "
243          << fsType;          << functionSpaceTypeAsString(fsType);
244      throw RipleyException(msg.str());      throw RipleyException(msg.str());
245  }  }
246    
247  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsCode, index_t id) const
248  {  {
249  #ifdef ESYS_MPI  #ifdef ESYS_MPI
250      if (fsCode == Nodes) {      if (fsCode != ReducedNodes) {
251          const index_t myFirst=getNumNodes()*m_mpiInfo->rank;          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
252          const index_t myLast=getNumNodes()*(m_mpiInfo->rank+1)-1;          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
253          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
254      } else      } else {
255          throw RipleyException("ownSample() only implemented for Nodes");          stringstream msg;
256            msg << "ownSample() not implemented for "
257                << functionSpaceTypeAsString(fsCode);
258            throw RipleyException(msg.str());
259        }
260  #else  #else
261      return true;      return true;
262  #endif  #endif
263  }  }
264    
265  void Rectangle::interpolateOnDomain(escript::Data& target,  void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const
                                     const escript::Data& in) const  
266  {  {
267      const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));      escript::Data& in = *const_cast<escript::Data*>(&cIn);
268      const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));      const dim_t numComp = in.getDataPointSize();
269      if (inDomain != *this)      const double h0 = m_l0/m_gNE0;
270          throw RipleyException("Illegal domain of interpolant");      const double h1 = m_l1/m_gNE1;
271      if (targetDomain != *this)      const double cx0 = -1./h0;
272          throw RipleyException("Illegal domain of interpolation target");      const double cx1 = -.78867513459481288225/h0;
273        const double cx2 = -.5/h0;
274      stringstream msg;      const double cx3 = -.21132486540518711775/h0;
275      msg << "interpolateOnDomain() not implemented for function space "      const double cx4 = .21132486540518711775/h0;
276          << in.getFunctionSpace().getTypeCode() << " -> "      const double cx5 = .5/h0;
277          << target.getFunctionSpace().getTypeCode();      const double cx6 = .78867513459481288225/h0;
278        const double cx7 = 1./h0;
279        const double cy0 = -1./h1;
280        const double cy1 = -.78867513459481288225/h1;
281        const double cy2 = -.5/h1;
282        const double cy3 = -.21132486540518711775/h1;
283        const double cy4 = .21132486540518711775/h1;
284        const double cy5 = .5/h1;
285        const double cy6 = .78867513459481288225/h1;
286        const double cy7 = 1./h1;
287    
288      switch (in.getFunctionSpace().getTypeCode()) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
289          case Nodes:          /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */
         case DegreesOfFreedom:  
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case Elements:  
                 {  
                     const double tmp0_2 = 0.62200846792814621559;  
                     const double tmp0_1 = 0.044658198738520451079;  
                     const double tmp0_0 = 0.16666666666666666667;  
                     const dim_t numComp = in.getDataPointSize();  
                     escript::Data* inn=const_cast<escript::Data*>(&in);  
290  #pragma omp parallel for  #pragma omp parallel for
291                      for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
292                          for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
293                              const register double* f_10 = inn->getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
294                              const register double* f_11 = inn->getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
295                              const register double* f_01 = inn->getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
296                              const register double* f_00 = inn->getSampleDataRO(INDEX2(k0,k1, m_N0));                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
297                              double* o = target.getSampleDataRW(INDEX2(k0,k1,m_NE0));                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
298                              for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
299                                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
300                                  o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
301                                  o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
302                                  o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
303                              } // close component loop i                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
304                          } // close k0 loop                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
305                      } // close k1 loop                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
306                      break;                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
307                    } /* end of component loop i */
308                } /* end of k0 loop */
309            } /* end of k1 loop */
310            /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
311        } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
312            /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */
313    #pragma omp parallel for
314            for (index_t k1=0; k1 < m_NE1; ++k1) {
315                for (index_t k0=0; k0 < m_NE0; ++k0) {
316                    const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
317                    const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
318                    const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
319                    const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
320                    double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
321                    for (index_t i=0; i < numComp; ++i) {
322                        o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
323                        o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
324                    } /* end of component loop i */
325                } /* end of k0 loop */
326            } /* end of k1 loop */
327            /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
328        } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
329    #pragma omp parallel
330            {
331                /*** GENERATOR SNIP_GRAD_FACES TOP */
332                if (m_faceOffset[0] > -1) {
333    #pragma omp for nowait
334                    for (index_t k1=0; k1 < m_NE1; ++k1) {
335                        const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
336                        const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
337                        const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
338                        const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
339                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
340                        for (index_t i=0; i < numComp; ++i) {
341                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
342                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
343                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
344                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
345                        } /* end of component loop i */
346                    } /* end of k1 loop */
347                } /* end of face 0 */
348                if (m_faceOffset[1] > -1) {
349    #pragma omp for nowait
350                    for (index_t k1=0; k1 < m_NE1; ++k1) {
351                        const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
352                        const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
353                        const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
354                        const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
355                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
356                        for (index_t i=0; i < numComp; ++i) {
357                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
358                            o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
359                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
360                            o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
361                        } /* end of component loop i */
362                    } /* end of k1 loop */
363                } /* end of face 1 */
364                if (m_faceOffset[2] > -1) {
365    #pragma omp for nowait
366                    for (index_t k0=0; k0 < m_NE0; ++k0) {
367                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
368                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
369                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
370                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
371                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
372                        for (index_t i=0; i < numComp; ++i) {
373                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
374                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
375                            o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
376                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
377                        } /* end of component loop i */
378                    } /* end of k0 loop */
379                } /* end of face 2 */
380                if (m_faceOffset[3] > -1) {
381    #pragma omp for nowait
382                    for (index_t k0=0; k0 < m_NE0; ++k0) {
383                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
384                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
385                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
386                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
387                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
388                        for (index_t i=0; i < numComp; ++i) {
389                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
390                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
391                            o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
392                            o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
393                        } /* end of component loop i */
394                    } /* end of k0 loop */
395                } /* end of face 3 */
396                /* GENERATOR SNIP_GRAD_FACES BOTTOM */
397            } // end of parallel section
398        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
399    #pragma omp parallel
400            {
401                /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
402                if (m_faceOffset[0] > -1) {
403    #pragma omp for nowait
404                    for (index_t k1=0; k1 < m_NE1; ++k1) {
405                        const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
406                        const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
407                        const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
408                        const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
409                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
410                        for (index_t i=0; i < numComp; ++i) {
411                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
412                            o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
413                        } /* end of component loop i */
414                    } /* end of k1 loop */
415                } /* end of face 0 */
416                if (m_faceOffset[1] > -1) {
417    #pragma omp for nowait
418                    for (index_t k1=0; k1 < m_NE1; ++k1) {
419                        const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
420                        const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
421                        const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
422                        const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
423                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
424                        for (index_t i=0; i < numComp; ++i) {
425                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
426                            o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
427                        } /* end of component loop i */
428                    } /* end of k1 loop */
429                } /* end of face 1 */
430                if (m_faceOffset[2] > -1) {
431    #pragma omp for nowait
432                    for (index_t k0=0; k0 < m_NE0; ++k0) {
433                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
434                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
435                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
436                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
437                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
438                        for (index_t i=0; i < numComp; ++i) {
439                            o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
440                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
441                        } /* end of component loop i */
442                    } /* end of k0 loop */
443                } /* end of face 2 */
444                if (m_faceOffset[3] > -1) {
445    #pragma omp for nowait
446                    for (index_t k0=0; k0 < m_NE0; ++k0) {
447                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
448                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
449                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
450                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
451                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
452                        for (index_t i=0; i < numComp; ++i) {
453                            o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
454                            o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
455                        } /* end of component loop i */
456                    } /* end of k0 loop */
457                } /* end of face 3 */
458                /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */
459            } // end of parallel section
460        } else {
461            stringstream msg;
462            msg << "setToGradient() not implemented for "
463                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
464            throw RipleyException(msg.str());
465        }
466    }
467    
468    void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
469    {
470        escript::Data& in = *const_cast<escript::Data*>(&arg);
471        const dim_t numComp = in.getDataPointSize();
472        const double h0 = m_l0/m_gNE0;
473        const double h1 = m_l1/m_gNE1;
474        if (arg.getFunctionSpace().getTypeCode() == Elements) {
475            const double w_0 = h0*h1/4.;
476    #pragma omp parallel
477            {
478                vector<double> int_local(numComp, 0);
479    #pragma omp for nowait
480                for (index_t k1 = 0; k1 < m_NE1; ++k1) {
481                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
482                        const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
483                        for (index_t i=0; i < numComp; ++i) {
484                            const register double f_0 = f[INDEX2(i,0,numComp)];
485                            const register double f_1 = f[INDEX2(i,1,numComp)];
486                            const register double f_2 = f[INDEX2(i,2,numComp)];
487                            const register double f_3 = f[INDEX2(i,3,numComp)];
488                            int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
489                        }  /* end of component loop i */
490                    } /* end of k0 loop */
491                } /* end of k1 loop */
492    
493    #pragma omp critical
494                for (index_t i=0; i<numComp; i++)
495                    integrals[i]+=int_local[i];
496            } // end of parallel section
497        } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
498            const double w_0 = h0*h1;
499    #pragma omp parallel
500            {
501                vector<double> int_local(numComp, 0);
502    #pragma omp for nowait
503                for (index_t k1 = 0; k1 < m_NE1; ++k1) {
504                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
505                        const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
506                        for (index_t i=0; i < numComp; ++i) {
507                            int_local[i]+=f[i]*w_0;
508                        }  /* end of component loop i */
509                    } /* end of k0 loop */
510                } /* end of k1 loop */
511    
512    #pragma omp critical
513                for (index_t i=0; i<numComp; i++)
514                    integrals[i]+=int_local[i];
515            } // end of parallel section
516        } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
517            const double w_0 = h0/2.;
518            const double w_1 = h1/2.;
519    #pragma omp parallel
520            {
521                vector<double> int_local(numComp, 0);
522                if (m_faceOffset[0] > -1) {
523    #pragma omp for nowait
524                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
525                        const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
526                        for (index_t i=0; i < numComp; ++i) {
527                            const register double f_0 = f[INDEX2(i,0,numComp)];
528                            const register double f_1 = f[INDEX2(i,1,numComp)];
529                            int_local[i]+=(f_0+f_1)*w_1;
530                        }  /* end of component loop i */
531                    } /* end of k1 loop */
532                }
533    
534                if (m_faceOffset[1] > -1) {
535    #pragma omp for nowait
536                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
537                        const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
538                        for (index_t i=0; i < numComp; ++i) {
539                            const register double f_0 = f[INDEX2(i,0,numComp)];
540                            const register double f_1 = f[INDEX2(i,1,numComp)];
541                            int_local[i]+=(f_0+f_1)*w_1;
542                        }  /* end of component loop i */
543                    } /* end of k1 loop */
544                }
545    
546                if (m_faceOffset[2] > -1) {
547    #pragma omp for nowait
548                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
549                        const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
550                        for (index_t i=0; i < numComp; ++i) {
551                            const register double f_0 = f[INDEX2(i,0,numComp)];
552                            const register double f_1 = f[INDEX2(i,1,numComp)];
553                            int_local[i]+=(f_0+f_1)*w_0;
554                        }  /* end of component loop i */
555                    } /* end of k0 loop */
556                }
557    
558                if (m_faceOffset[3] > -1) {
559    #pragma omp for nowait
560                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
561                        const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
562                        for (index_t i=0; i < numComp; ++i) {
563                            const register double f_0 = f[INDEX2(i,0,numComp)];
564                            const register double f_1 = f[INDEX2(i,1,numComp)];
565                            int_local[i]+=(f_0+f_1)*w_0;
566                        }  /* end of component loop i */
567                    } /* end of k0 loop */
568                }
569    
570    #pragma omp critical
571                for (index_t i=0; i<numComp; i++)
572                    integrals[i]+=int_local[i];
573            } // end of parallel section
574        } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
575    #pragma omp parallel
576            {
577                vector<double> int_local(numComp, 0);
578                if (m_faceOffset[0] > -1) {
579    #pragma omp for nowait
580                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
581                        const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
582                        for (index_t i=0; i < numComp; ++i) {
583                            int_local[i]+=f[i]*h1;
584                        }  /* end of component loop i */
585                    } /* end of k1 loop */
586                }
587    
588                if (m_faceOffset[1] > -1) {
589    #pragma omp for nowait
590                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
591                        const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
592                        for (index_t i=0; i < numComp; ++i) {
593                            int_local[i]+=f[i]*h1;
594                        }  /* end of component loop i */
595                    } /* end of k1 loop */
596                }
597    
598                if (m_faceOffset[2] > -1) {
599    #pragma omp for nowait
600                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
601                        const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
602                        for (index_t i=0; i < numComp; ++i) {
603                            int_local[i]+=f[i]*h0;
604                        }  /* end of component loop i */
605                    } /* end of k0 loop */
606                }
607    
608                if (m_faceOffset[3] > -1) {
609    #pragma omp for nowait
610                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
611                        const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
612                        for (index_t i=0; i < numComp; ++i) {
613                            int_local[i]+=f[i]*h0;
614                        }  /* end of component loop i */
615                    } /* end of k0 loop */
616                }
617    
618    #pragma omp critical
619                for (index_t i=0; i<numComp; i++)
620                    integrals[i]+=int_local[i];
621            } // end of parallel section
622        } else {
623            stringstream msg;
624            msg << "setToIntegrals() not implemented for "
625                << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
626            throw RipleyException(msg.str());
627        }
628    }
629    
630    void Rectangle::setToNormal(escript::Data& out) const
631    {
632        if (out.getFunctionSpace().getTypeCode() == FaceElements) {
633    #pragma omp parallel
634            {
635                if (m_faceOffset[0] > -1) {
636    #pragma omp for nowait
637                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
638                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
639                        // set vector at two quadrature points
640                        *o++ = -1.;
641                        *o++ = 0.;
642                        *o++ = -1.;
643                        *o = 0.;
644                  }                  }
645                }
646    
647                  case DegreesOfFreedom:              if (m_faceOffset[1] > -1) {
648                      target=in;  #pragma omp for nowait
649                      break;                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
650                                            double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
651                  default:                      // set vector at two quadrature points
652                      throw RipleyException(msg.str());                      *o++ = 1.;
653                        *o++ = 0.;
654                        *o++ = 1.;
655                        *o = 0.;
656                    }
657              }              }
658          default:  
659              throw RipleyException(msg.str());              if (m_faceOffset[2] > -1) {
660    #pragma omp for nowait
661                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
662                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
663                        // set vector at two quadrature points
664                        *o++ = 0.;
665                        *o++ = -1.;
666                        *o++ = 0.;
667                        *o = -1.;
668                    }
669                }
670    
671                if (m_faceOffset[3] > -1) {
672    #pragma omp for nowait
673                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
674                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
675                        // set vector at two quadrature points
676                        *o++ = 0.;
677                        *o++ = 1.;
678                        *o++ = 0.;
679                        *o = 1.;
680                    }
681                }
682            } // end of parallel section
683        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
684    #pragma omp parallel
685            {
686                if (m_faceOffset[0] > -1) {
687    #pragma omp for nowait
688                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
689                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
690                        *o++ = -1.;
691                        *o = 0.;
692                    }
693                }
694    
695                if (m_faceOffset[1] > -1) {
696    #pragma omp for nowait
697                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
698                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
699                        *o++ = 1.;
700                        *o = 0.;
701                    }
702                }
703    
704                if (m_faceOffset[2] > -1) {
705    #pragma omp for nowait
706                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
707                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
708                        *o++ = 0.;
709                        *o = -1.;
710                    }
711                }
712    
713                if (m_faceOffset[3] > -1) {
714    #pragma omp for nowait
715                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
716                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
717                        *o++ = 0.;
718                        *o = 1.;
719                    }
720                }
721            } // end of parallel section
722    
723        } else {
724            stringstream msg;
725            msg << "setToNormal() not implemented for "
726                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
727            throw RipleyException(msg.str());
728      }      }
729  }  }
730    
# Line 313  Paso_SystemMatrixPattern* Rectangle::get Line 734  Paso_SystemMatrixPattern* Rectangle::get
734      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
735          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
736    
737  /*      // connector
738      // create distribution      RankVector neighbour;
739      IndexVector dist;      IndexVector offsetInShared(1,0);
740      for (index_t i=0; i<m_mpiInfo->size+1; i++)      IndexVector sendShared, recvShared;
741          dist.push_back(i*getNumNodes());      const IndexVector faces=getNumFacesPerBoundary();
742      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      const index_t left = (m_offset0==0 ? 0 : 1);
743              &dist[0], 1, 0);      const index_t bottom = (m_offset1==0 ? 0 : 1);
744        // corner node from bottom-left
745      // connectors      if (faces[0] == 0 && faces[2] == 0) {
746      dim_t numNeighbours = 0;          neighbour.push_back(m_mpiInfo->rank-m_NX-1);
747      RankVector neighbour(numNeighbours);          offsetInShared.push_back(offsetInShared.back()+1);
748      IndexVector offsetInShared(numNeighbours+1);          sendShared.push_back(m_nodeId[m_N0+left]);
749      IndexVector shared(offsetInShared[numNeighbours]);          recvShared.push_back(m_nodeId[0]);
750        }
751        // bottom edge
752        if (faces[2] == 0) {
753            neighbour.push_back(m_mpiInfo->rank-m_NX);
754            offsetInShared.push_back(offsetInShared.back()+m_N0-left);
755            for (dim_t i=left; i<m_N0; i++) {
756                // easy case, we know the neighbour id's
757                sendShared.push_back(m_nodeId[m_N0+i]);
758                recvShared.push_back(m_nodeId[i]);
759            }
760        }
761        // corner node from bottom-right
762        if (faces[1] == 0 && faces[2] == 0) {
763            neighbour.push_back(m_mpiInfo->rank-m_NX+1);
764            const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);
765            const index_t N1=(neighbour.back()/m_NX == 0 ? m_N1 : m_N1-1);
766            const index_t first=m_nodeDistribution[neighbour.back()];
767            offsetInShared.push_back(offsetInShared.back()+1);
768            sendShared.push_back(m_nodeId[(bottom+1)*m_N0-1]);
769            recvShared.push_back(first+N0*(N1-1));
770        }
771        // left edge
772        if (faces[0] == 0) {
773            neighbour.push_back(m_mpiInfo->rank-1);
774            offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);
775            for (dim_t i=bottom; i<m_N1; i++) {
776                // easy case, we know the neighbour id's
777                sendShared.push_back(m_nodeId[i*m_N0+1]);
778                recvShared.push_back(m_nodeId[i*m_N0]);
779            }
780        }
781        // right edge
782        if (faces[1] == 0) {
783            neighbour.push_back(m_mpiInfo->rank+1);
784            const index_t rightN0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);
785            const index_t first=m_nodeDistribution[neighbour.back()];
786            offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);
787            for (dim_t i=bottom; i<m_N1; i++) {
788                sendShared.push_back(m_nodeId[(i+1)*m_N0-1]);
789                recvShared.push_back(first+rightN0*(i-bottom));
790            }
791        }
792        // corner node from top-left
793        if (faces[0] == 0 && faces[3] == 0) {
794            neighbour.push_back(m_mpiInfo->rank+m_NX-1);
795            const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);
796            const index_t first=m_nodeDistribution[neighbour.back()];
797            offsetInShared.push_back(offsetInShared.back()+1);
798            sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+left]);
799            recvShared.push_back(first+N0-1);
800        }
801        // top edge
802        if (faces[3] == 0) {
803            neighbour.push_back(m_mpiInfo->rank+m_NX);
804            const index_t first=m_nodeDistribution[neighbour.back()];
805            offsetInShared.push_back(offsetInShared.back()+m_N0-left);
806            for (dim_t i=left; i<m_N0; i++) {
807                sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+i]);
808                recvShared.push_back(first+i-left);
809            }
810        }
811        // corner node from top-right
812        if (faces[1] == 0 && faces[3] == 0) {
813            neighbour.push_back(m_mpiInfo->rank+m_NX+1);
814            const index_t first=m_nodeDistribution[neighbour.back()];
815            offsetInShared.push_back(offsetInShared.back()+1);
816            sendShared.push_back(m_nodeId[m_N0*m_N1-1]);
817            recvShared.push_back(first);
818        }
819        const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];
820        /*
821        cout << "--- rcv_shcomp ---" << endl;
822        cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
823        for (size_t i=0; i<neighbour.size(); i++) {
824            cout << "neighbor[" << i << "]=" << neighbour[i]
825                << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
826        }
827        for (size_t i=0; i<recvShared.size(); i++) {
828            cout << "shared[" << i << "]=" << recvShared[i] << endl;
829        }
830        cout << "--- snd_shcomp ---" << endl;
831        for (size_t i=0; i<sendShared.size(); i++) {
832            cout << "shared[" << i << "]=" << sendShared[i] << endl;
833        }
834        */
835    
836      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
837              getNumNodes(), numNeighbours, &neighbour[0], &shared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
838              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo);
839      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
840              getNumNodes(), numNeighbours, &neighbour[0], &shared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
841              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo);
842      Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);      Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
843        Paso_SharedComponents_free(snd_shcomp);
844        Paso_SharedComponents_free(rcv_shcomp);
845    
846      // create patterns      // create patterns
847      dim_t M=0, N=0;      dim_t M, N;
848      int* ptr=NULL;      IndexVector ptr(1,0);
849      index_t* index=NULL;      IndexVector index;
850    
851        // main pattern
852        for (index_t i=0; i<numDOF; i++) {
853            // always add the node itself
854            index.push_back(i);
855            const int num=insertNeighbours(index, i);
856            ptr.push_back(ptr.back()+num+1);
857        }
858        M=N=ptr.size()-1;
859        // paso will manage the memory
860        index_t* indexC = MEMALLOC(index.size(),index_t);
861        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
862        copy(index.begin(), index.end(), indexC);
863        copy(ptr.begin(), ptr.end(), ptrC);
864      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
865              M, N, ptr, index);              M, N, ptrC, indexC);
866      Paso_Pattern *colCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
867              M, N, ptr, index);      /*
868      Paso_Pattern *rowCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      cout << "--- main_pattern ---" << endl;
869              M, N, ptr, index);      cout << "M=" << M << ", N=" << N << endl;
870        for (size_t i=0; i<ptr.size(); i++) {
871            cout << "ptr[" << i << "]=" << ptr[i] << endl;
872        }
873        for (size_t i=0; i<index.size(); i++) {
874            cout << "index[" << i << "]=" << index[i] << endl;
875        }
876        */
877    
878        ptr.clear();
879        index.clear();
880    
881        // column & row couple patterns
882        Paso_Pattern *colCouplePattern, *rowCouplePattern;
883        generateCouplePatterns(&colCouplePattern, &rowCouplePattern);
884    
885        // allocate paso distribution
886        Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
887                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
888    
889      Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(      Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(
890              MATRIX_FORMAT_DEFAULT, distribution, distribution,              MATRIX_FORMAT_DEFAULT, distribution, distribution,
# Line 354  Paso_SystemMatrixPattern* Rectangle::get Line 894  Paso_SystemMatrixPattern* Rectangle::get
894      Paso_Pattern_free(colCouplePattern);      Paso_Pattern_free(colCouplePattern);
895      Paso_Pattern_free(rowCouplePattern);      Paso_Pattern_free(rowCouplePattern);
896      Paso_Distribution_free(distribution);      Paso_Distribution_free(distribution);
     Paso_SharedComponents_free(snd_shcomp);  
     Paso_SharedComponents_free(rcv_shcomp);  
897      return pattern;      return pattern;
 */  
     throw RipleyException("getPattern() not implemented");  
898  }  }
899    
900  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
# Line 425  pair<double,double> Rectangle::getFirstC Line 961  pair<double,double> Rectangle::getFirstC
961  //protected  //protected
962  dim_t Rectangle::getNumFaceElements() const  dim_t Rectangle::getNumFaceElements() const
963  {  {
964        const IndexVector faces = getNumFacesPerBoundary();
965      dim_t n=0;      dim_t n=0;
966      //left      for (size_t i=0; i<faces.size(); i++)
967      if (m_offset0==0)          n+=faces[i];
         n+=m_NE1;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         n+=m_NE1;  
     //bottom  
     if (m_offset1==0)  
         n+=m_NE0;  
     //top  
     if (m_mpiInfo->rank/m_NX==m_NY-1)  
         n+=m_NE0;  
   
968      return n;      return n;
969  }  }
970    
# Line 518  void Rectangle::populateSampleIds() Line 1044  void Rectangle::populateSampleIds()
1044      }      }
1045      if (left>0 && bottom>0) {      if (left>0 && bottom>0) {
1046          const int neighbour=m_mpiInfo->rank-m_NX-1;          const int neighbour=m_mpiInfo->rank-m_NX-1;
1047          const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);          const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1048          const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);          const index_t N1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);
1049          m_nodeId[0]=m_nodeDistribution[neighbour]+bottomN1*bottomN0-1;          m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1-1;
1050      }      }
1051    
1052      // the rest of the id's are contiguous      // the rest of the id's are contiguous
# Line 531  void Rectangle::populateSampleIds() Line 1057  void Rectangle::populateSampleIds()
1057              m_nodeId[i0+i1*m_N0] = firstId+i0-left+(i1-bottom)*(m_N0-left);              m_nodeId[i0+i1*m_N0] = firstId+i0-left+(i1-bottom)*(m_N0-left);
1058          }          }
1059      }      }
1060        m_nodeTags.assign(getNumNodes(), 0);
1061      for (dim_t k=0; k<m_nodeDistribution.size(); k++) {      updateTagsInUse(Nodes);
         cout << m_nodeDistribution[k] << endl;  
     }  
     cout<<endl;  
     for (dim_t k=0; k<getNumNodes(); k++) {  
         cout << m_nodeId[k] << endl;  
     }  
1062    
1063      // elements      // elements
1064      m_elementId.resize(getNumElements());      m_elementId.resize(getNumElements());
# Line 546  void Rectangle::populateSampleIds() Line 1066  void Rectangle::populateSampleIds()
1066      for (dim_t k=0; k<getNumElements(); k++) {      for (dim_t k=0; k<getNumElements(); k++) {
1067          m_elementId[k]=k;          m_elementId[k]=k;
1068      }      }
1069        m_elementTags.assign(getNumElements(), 0);
1070        updateTagsInUse(Elements);
1071    
1072      // face elements      // face element id's
1073      m_faceId.resize(getNumFaceElements());      m_faceId.resize(getNumFaceElements());
1074  #pragma omp parallel for  #pragma omp parallel for
1075      for (dim_t k=0; k<getNumFaceElements(); k++) {      for (dim_t k=0; k<getNumFaceElements(); k++) {
1076          m_faceId[k]=k;          m_faceId[k]=k;
1077      }      }
1078    
1079        // generate face offset vector and set face tags
1080        const IndexVector facesPerEdge = getNumFacesPerBoundary();
1081        const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
1082        const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
1083        m_faceOffset.assign(facesPerEdge.size(), -1);
1084        m_faceTags.clear();
1085        index_t offset=0;
1086        for (size_t i=0; i<facesPerEdge.size(); i++) {
1087            if (facesPerEdge[i]>0) {
1088                m_faceOffset[i]=offset;
1089                offset+=facesPerEdge[i];
1090                m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);
1091            }
1092        }
1093        setTagMap("left", LEFT);
1094        setTagMap("right", RIGHT);
1095        setTagMap("bottom", BOTTOM);
1096        setTagMap("top", TOP);
1097        updateTagsInUse(FaceElements);
1098    }
1099    
1100    //private
1101    int Rectangle::insertNeighbours(IndexVector& index, index_t node) const
1102    {
1103        const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);
1104        const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);
1105        const int x=node%myN0;
1106        const int y=node/myN0;
1107        int num=0;
1108        if (y>0) {
1109            if (x>0) {
1110                // bottom-left
1111                index.push_back(node-myN0-1);
1112                num++;
1113            }
1114            // bottom
1115            index.push_back(node-myN0);
1116            num++;
1117            if (x<myN0-1) {
1118                // bottom-right
1119                index.push_back(node-myN0+1);
1120                num++;
1121            }
1122        }
1123        if (x<myN0-1) {
1124            // right
1125            index.push_back(node+1);
1126            num++;
1127            if (y<myN1-1) {
1128                // top-right
1129                index.push_back(node+myN0+1);
1130                num++;
1131            }
1132        }
1133        if (y<myN1-1) {
1134            // top
1135            index.push_back(node+myN0);
1136            num++;
1137            if (x>0) {
1138                // top-left
1139                index.push_back(node+myN0-1);
1140                num++;
1141            }
1142        }
1143        if (x>0) {
1144            // left
1145            index.push_back(node-1);
1146            num++;
1147        }
1148    
1149        return num;
1150    }
1151    
1152    //private
1153    void Rectangle::generateCouplePatterns(Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
1154    {
1155        IndexVector ptr(1,0);
1156        IndexVector index;
1157        const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);
1158        const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);
1159        const IndexVector faces=getNumFacesPerBoundary();
1160    
1161        // bottom edge
1162        dim_t n=0;
1163        if (faces[0] == 0) {
1164            index.push_back(2*(myN0+myN1+1));
1165            index.push_back(2*(myN0+myN1+1)+1);
1166            n+=2;
1167            if (faces[2] == 0) {
1168                index.push_back(0);
1169                index.push_back(1);
1170                index.push_back(2);
1171                n+=3;
1172            }
1173        } else if (faces[2] == 0) {
1174            index.push_back(1);
1175            index.push_back(2);
1176            n+=2;
1177        }
1178        // n=neighbours of bottom-left corner node
1179        ptr.push_back(ptr.back()+n);
1180        n=0;
1181        if (faces[2] == 0) {
1182            for (dim_t i=1; i<myN0-1; i++) {
1183                index.push_back(i);
1184                index.push_back(i+1);
1185                index.push_back(i+2);
1186                ptr.push_back(ptr.back()+3);
1187            }
1188            index.push_back(myN0-1);
1189            index.push_back(myN0);
1190            n+=2;
1191            if (faces[1] == 0) {
1192                index.push_back(myN0+1);
1193                index.push_back(myN0+2);
1194                index.push_back(myN0+3);
1195                n+=3;
1196            }
1197        } else {
1198            for (dim_t i=1; i<myN0-1; i++) {
1199                ptr.push_back(ptr.back());
1200            }
1201            if (faces[1] == 0) {
1202                index.push_back(myN0+2);
1203                index.push_back(myN0+3);
1204                n+=2;
1205            }
1206        }
1207        // n=neighbours of bottom-right corner node
1208        ptr.push_back(ptr.back()+n);
1209    
1210        // 2nd row to 2nd last row
1211        for (dim_t i=1; i<myN1-1; i++) {
1212            // left edge
1213            if (faces[0] == 0) {
1214                index.push_back(2*(myN0+myN1+2)-i);
1215                index.push_back(2*(myN0+myN1+2)-i-1);
1216                index.push_back(2*(myN0+myN1+2)-i-2);
1217                ptr.push_back(ptr.back()+3);
1218            } else {
1219                ptr.push_back(ptr.back());
1220            }
1221            for (dim_t j=1; j<myN0-1; j++) {
1222                ptr.push_back(ptr.back());
1223            }
1224            // right edge
1225            if (faces[1] == 0) {
1226                index.push_back(myN0+i+1);
1227                index.push_back(myN0+i+2);
1228                index.push_back(myN0+i+3);
1229                ptr.push_back(ptr.back()+3);
1230            } else {
1231                ptr.push_back(ptr.back());
1232            }
1233        }
1234    
1235        // top edge
1236        n=0;
1237        if (faces[0] == 0) {
1238            index.push_back(2*myN0+myN1+5);
1239            index.push_back(2*myN0+myN1+4);
1240            n+=2;
1241            if (faces[3] == 0) {
1242                index.push_back(2*myN0+myN1+3);
1243                index.push_back(2*myN0+myN1+2);
1244                index.push_back(2*myN0+myN1+1);
1245                n+=3;
1246            }
1247        } else if (faces[3] == 0) {
1248            index.push_back(2*myN0+myN1+2);
1249            index.push_back(2*myN0+myN1+1);
1250            n+=2;
1251        }
1252        // n=neighbours of top-left corner node
1253        ptr.push_back(ptr.back()+n);
1254        n=0;
1255        if (faces[3] == 0) {
1256            for (dim_t i=1; i<myN0-1; i++) {
1257                index.push_back(2*myN0+myN1+i+1);
1258                index.push_back(2*myN0+myN1+i);
1259                index.push_back(2*myN0+myN1+i-1);
1260                ptr.push_back(ptr.back()+3);
1261            }
1262            index.push_back(myN0+myN1+4);
1263            index.push_back(myN0+myN1+3);
1264            n+=2;
1265            if (faces[1] == 0) {
1266                index.push_back(myN0+myN1+2);
1267                index.push_back(myN0+myN1+1);
1268                index.push_back(myN0+myN1);
1269                n+=3;
1270            }
1271        } else {
1272            for (dim_t i=1; i<myN0-1; i++) {
1273                ptr.push_back(ptr.back());
1274            }
1275            if (faces[1] == 0) {
1276                index.push_back(myN0+myN1+1);
1277                index.push_back(myN0+myN1);
1278                n+=2;
1279            }
1280        }
1281        // n=neighbours of top-right corner node
1282        ptr.push_back(ptr.back()+n);
1283    
1284        dim_t M=ptr.size()-1;
1285        map<index_t,index_t> idMap;
1286        dim_t N=0;
1287        for (IndexVector::iterator it=index.begin(); it!=index.end(); it++) {
1288            if (idMap.find(*it)==idMap.end()) {
1289                idMap[*it]=N;
1290                *it=N++;
1291            } else {
1292                *it=idMap[*it];
1293            }
1294        }
1295    
1296        /*
1297        cout << "--- colCouple_pattern ---" << endl;
1298        cout << "M=" << M << ", N=" << N << endl;
1299        for (size_t i=0; i<ptr.size(); i++) {
1300            cout << "ptr[" << i << "]=" << ptr[i] << endl;
1301        }
1302        for (size_t i=0; i<index.size(); i++) {
1303            cout << "index[" << i << "]=" << index[i] << endl;
1304        }
1305        */
1306    
1307        // now build the row couple pattern
1308        IndexVector ptr2(1,0);
1309        IndexVector index2;
1310        for (dim_t id=0; id<N; id++) {
1311            n=0;
1312            for (dim_t i=0; i<M; i++) {
1313                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
1314                    if (index[j]==id) {
1315                        index2.push_back(i);
1316                        n++;
1317                        break;
1318                    }
1319                }
1320            }
1321            ptr2.push_back(ptr2.back()+n);
1322        }
1323    
1324        /*
1325        cout << "--- rowCouple_pattern ---" << endl;
1326        cout << "M=" << N << ", N=" << M << endl;
1327        for (size_t i=0; i<ptr2.size(); i++) {
1328            cout << "ptr[" << i << "]=" << ptr2[i] << endl;
1329        }
1330        for (size_t i=0; i<index2.size(); i++) {
1331            cout << "index[" << i << "]=" << index2[i] << endl;
1332        }
1333        */
1334    
1335        // paso will manage the memory
1336        index_t* indexC = MEMALLOC(index.size(), index_t);
1337        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
1338        copy(index.begin(), index.end(), indexC);
1339        copy(ptr.begin(), ptr.end(), ptrC);
1340        *colPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
1341    
1342        // paso will manage the memory
1343        indexC = MEMALLOC(index2.size(), index_t);
1344        ptrC = MEMALLOC(ptr2.size(), index_t);
1345        copy(index2.begin(), index2.end(), indexC);
1346        copy(ptr2.begin(), ptr2.end(), ptrC);
1347        *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);
1348    }
1349    
1350    //protected
1351    void Rectangle::interpolateNodesOnElements(escript::Data& out,
1352                                            escript::Data& in, bool reduced) const
1353    {
1354        const dim_t numComp = in.getDataPointSize();
1355        if (reduced) {
1356            /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */
1357            const double c0 = .25;
1358    #pragma omp parallel for
1359            for (index_t k1=0; k1 < m_NE1; ++k1) {
1360                for (index_t k0=0; k0 < m_NE0; ++k0) {
1361                    const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1362                    const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1363                    const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1364                    const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1365                    double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1366                    for (index_t i=0; i < numComp; ++i) {
1367                        o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1368                    } /* end of component loop i */
1369                } /* end of k0 loop */
1370            } /* end of k1 loop */
1371            /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */
1372        } else {
1373            /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
1374            const double c0 = .16666666666666666667;
1375            const double c1 = .044658198738520451079;
1376            const double c2 = .62200846792814621559;
1377    #pragma omp parallel for
1378            for (index_t k1=0; k1 < m_NE1; ++k1) {
1379                for (index_t k0=0; k0 < m_NE0; ++k0) {
1380                    const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
1381                    const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
1382                    const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
1383                    const register double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));
1384                    double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1385                    for (index_t i=0; i < numComp; ++i) {
1386                        o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);
1387                        o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);
1388                        o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);
1389                        o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);
1390                    } /* end of component loop i */
1391                } /* end of k0 loop */
1392            } /* end of k1 loop */
1393            /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
1394        }
1395    }
1396    
1397    //protected
1398    void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,
1399                                            bool reduced) const
1400    {
1401        const dim_t numComp = in.getDataPointSize();
1402        if (reduced) {
1403            const double c0 = .5;
1404    #pragma omp parallel
1405            {
1406                /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */
1407                if (m_faceOffset[0] > -1) {
1408    #pragma omp for nowait
1409                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1410                        const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1411                        const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1412                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1413                        for (index_t i=0; i < numComp; ++i) {
1414                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
1415                        } /* end of component loop i */
1416                    } /* end of k1 loop */
1417                } /* end of face 0 */
1418                if (m_faceOffset[1] > -1) {
1419    #pragma omp for nowait
1420                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1421                        const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1422                        const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1423                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1424                        for (index_t i=0; i < numComp; ++i) {
1425                            o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
1426                        } /* end of component loop i */
1427                    } /* end of k1 loop */
1428                } /* end of face 1 */
1429                if (m_faceOffset[2] > -1) {
1430    #pragma omp for nowait
1431                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1432                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1433                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1434                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1435                        for (index_t i=0; i < numComp; ++i) {
1436                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
1437                        } /* end of component loop i */
1438                    } /* end of k0 loop */
1439                } /* end of face 2 */
1440                if (m_faceOffset[3] > -1) {
1441    #pragma omp for nowait
1442                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1443                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1444                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1445                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1446                        for (index_t i=0; i < numComp; ++i) {
1447                            o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1448                        } /* end of component loop i */
1449                    } /* end of k0 loop */
1450                } /* end of face 3 */
1451                /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */
1452            } // end of parallel section
1453        } else {
1454            const double c0 = 0.21132486540518711775;
1455            const double c1 = 0.78867513459481288225;
1456    #pragma omp parallel
1457            {
1458                /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */
1459                if (m_faceOffset[0] > -1) {
1460    #pragma omp for nowait
1461                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1462                        const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1463                        const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1464                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1465                        for (index_t i=0; i < numComp; ++i) {
1466                            o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;
1467                            o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_01[i]*c1;
1468                        } /* end of component loop i */
1469                    } /* end of k1 loop */
1470                } /* end of face 0 */
1471                if (m_faceOffset[1] > -1) {
1472    #pragma omp for nowait
1473                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1474                        const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1475                        const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1476                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1477                        for (index_t i=0; i < numComp; ++i) {
1478                            o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;
1479                            o[INDEX2(i,numComp,1)] = f_10[i]*c0 + f_11[i]*c1;
1480                        } /* end of component loop i */
1481                    } /* end of k1 loop */
1482                } /* end of face 1 */
1483                if (m_faceOffset[2] > -1) {
1484    #pragma omp for nowait
1485                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1486                        const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1487                        const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1488                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1489                        for (index_t i=0; i < numComp; ++i) {
1490                            o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;
1491                            o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_10[i]*c1;
1492                        } /* end of component loop i */
1493                    } /* end of k0 loop */
1494                } /* end of face 2 */
1495                if (m_faceOffset[3] > -1) {
1496    #pragma omp for nowait
1497                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1498                        const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1499                        const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1500                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1501                        for (index_t i=0; i < numComp; ++i) {
1502                            o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;
1503                            o[INDEX2(i,numComp,1)] = f_01[i]*c0 + f_11[i]*c1;
1504                        } /* end of component loop i */
1505                    } /* end of k0 loop */
1506                } /* end of face 3 */
1507                /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
1508            } // end of parallel section
1509        }
1510    }
1511    
1512    //protected
1513    void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1514            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1515            const escript::Data& C, const escript::Data& D,
1516            const escript::Data& X, const escript::Data& Y,
1517            const escript::Data& d, const escript::Data& y) const
1518    {
1519        const double h0 = m_l0/m_gNE0;
1520        const double h1 = m_l1/m_gNE1;
1521        /* GENERATOR SNIP_PDE_SINGLE_PRE TOP */
1522        const double w0 = -0.1555021169820365539*h1/h0;
1523        const double w1 = 0.041666666666666666667;
1524        const double w10 = -0.041666666666666666667*h0/h1;
1525        const double w11 = 0.1555021169820365539*h1/h0;
1526        const double w12 = 0.1555021169820365539*h0/h1;
1527        const double w13 = 0.01116454968463011277*h0/h1;
1528        const double w14 = 0.01116454968463011277*h1/h0;
1529        const double w15 = 0.041666666666666666667*h1/h0;
1530        const double w16 = -0.01116454968463011277*h0/h1;
1531        const double w17 = -0.1555021169820365539*h0/h1;
1532        const double w18 = -0.33333333333333333333*h1/h0;
1533        const double w19 = 0.25000000000000000000;
1534        const double w2 = -0.15550211698203655390;
1535        const double w20 = -0.25000000000000000000;
1536        const double w21 = 0.16666666666666666667*h0/h1;
1537        const double w22 = -0.16666666666666666667*h1/h0;
1538        const double w23 = -0.16666666666666666667*h0/h1;
1539        const double w24 = 0.33333333333333333333*h1/h0;
1540        const double w25 = 0.33333333333333333333*h0/h1;
1541        const double w26 = 0.16666666666666666667*h1/h0;
1542        const double w27 = -0.33333333333333333333*h0/h1;
1543        const double w28 = -0.032861463941450536761*h1;
1544        const double w29 = -0.032861463941450536761*h0;
1545        const double w3 = 0.041666666666666666667*h0/h1;
1546        const double w30 = -0.12264065304058601714*h1;
1547        const double w31 = -0.0023593469594139828636*h1;
1548        const double w32 = -0.008805202725216129906*h0;
1549        const double w33 = -0.008805202725216129906*h1;
1550        const double w34 = 0.032861463941450536761*h1;
1551        const double w35 = 0.008805202725216129906*h1;
1552        const double w36 = 0.008805202725216129906*h0;
1553        const double w37 = 0.0023593469594139828636*h1;
1554        const double w38 = 0.12264065304058601714*h1;
1555        const double w39 = 0.032861463941450536761*h0;
1556        const double w4 = 0.15550211698203655390;
1557        const double w40 = -0.12264065304058601714*h0;
1558        const double w41 = -0.0023593469594139828636*h0;
1559        const double w42 = 0.0023593469594139828636*h0;
1560        const double w43 = 0.12264065304058601714*h0;
1561        const double w44 = -0.16666666666666666667*h1;
1562        const double w45 = -0.083333333333333333333*h0;
1563        const double w46 = 0.083333333333333333333*h1;
1564        const double w47 = 0.16666666666666666667*h1;
1565        const double w48 = 0.083333333333333333333*h0;
1566        const double w49 = -0.16666666666666666667*h0;
1567        const double w5 = -0.041666666666666666667;
1568        const double w50 = 0.16666666666666666667*h0;
1569        const double w51 = -0.083333333333333333333*h1;
1570        const double w52 = 0.025917019497006092316*h0*h1;
1571        const double w53 = 0.0018607582807716854616*h0*h1;
1572        const double w54 = 0.0069444444444444444444*h0*h1;
1573        const double w55 = 0.09672363354357992482*h0*h1;
1574        const double w56 = 0.00049858867864229740201*h0*h1;
1575        const double w57 = 0.055555555555555555556*h0*h1;
1576        const double w58 = 0.027777777777777777778*h0*h1;
1577        const double w59 = 0.11111111111111111111*h0*h1;
1578        const double w6 = -0.01116454968463011277*h1/h0;
1579        const double w60 = -0.19716878364870322056*h1;
1580        const double w61 = -0.19716878364870322056*h0;
1581        const double w62 = -0.052831216351296779436*h0;
1582        const double w63 = -0.052831216351296779436*h1;
1583        const double w64 = 0.19716878364870322056*h1;
1584        const double w65 = 0.052831216351296779436*h1;
1585        const double w66 = 0.19716878364870322056*h0;
1586        const double w67 = 0.052831216351296779436*h0;
1587        const double w68 = -0.5*h1;
1588        const double w69 = -0.5*h0;
1589        const double w7 = 0.011164549684630112770;
1590        const double w70 = 0.5*h1;
1591        const double w71 = 0.5*h0;
1592        const double w72 = 0.1555021169820365539*h0*h1;
1593        const double w73 = 0.041666666666666666667*h0*h1;
1594        const double w74 = 0.01116454968463011277*h0*h1;
1595        const double w75 = 0.25*h0*h1;
1596        const double w8 = -0.011164549684630112770;
1597        const double w9 = -0.041666666666666666667*h1/h0;
1598        /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */
1599    
1600        rhs.requireWrite();
1601    #pragma omp parallel
1602        {
1603            for (index_t k1_0=0; k1_0<2; k1_0++) { // coloring
1604    #pragma omp for
1605                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
1606                    for (index_t k0=0; k0<m_NE0; ++k0)  {
1607                        bool add_EM_S=false;
1608                        bool add_EM_F=false;
1609                        vector<double> EM_S(4*4, 0);
1610                        vector<double> EM_F(4, 0);
1611                        const index_t e = k0 + m_NE0*k1;
1612                        /* GENERATOR SNIP_PDE_SINGLE TOP */
1613                        ///////////////
1614                        // process A //
1615                        ///////////////
1616                        if (!A.isEmpty()) {
1617                            add_EM_S=true;
1618                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1619                            if (A.actsExpanded()) {
1620                                const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
1621                                const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1622                                const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1623                                const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1624                                const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
1625                                const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1626                                const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1627                                const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1628                                const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
1629                                const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1630                                const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1631                                const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1632                                const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
1633                                const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1634                                const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1635                                const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1636                                const register double tmp4_0 = A_10_1 + A_10_2;
1637                                const register double tmp12_0 = A_11_0 + A_11_2;
1638                                const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1639                                const register double tmp10_0 = A_01_3 + A_10_3;
1640                                const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
1641                                const register double tmp0_0 = A_01_0 + A_01_3;
1642                                const register double tmp13_0 = A_01_2 + A_10_1;
1643                                const register double tmp3_0 = A_00_2 + A_00_3;
1644                                const register double tmp11_0 = A_11_1 + A_11_3;
1645                                const register double tmp18_0 = A_01_1 + A_10_1;
1646                                const register double tmp1_0 = A_00_0 + A_00_1;
1647                                const register double tmp15_0 = A_01_1 + A_10_2;
1648                                const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1649                                const register double tmp16_0 = A_10_0 + A_10_3;
1650                                const register double tmp6_0 = A_01_3 + A_10_0;
1651                                const register double tmp17_0 = A_01_1 + A_01_2;
1652                                const register double tmp9_0 = A_01_0 + A_10_0;
1653                                const register double tmp7_0 = A_01_0 + A_10_3;
1654                                const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1655                                const register double tmp19_0 = A_01_2 + A_10_2;
1656                                const register double tmp14_1 = A_10_0*w8;
1657                                const register double tmp23_1 = tmp3_0*w14;
1658                                const register double tmp35_1 = A_01_0*w8;
1659                                const register double tmp54_1 = tmp13_0*w8;
1660                                const register double tmp20_1 = tmp9_0*w4;
1661                                const register double tmp25_1 = tmp12_0*w12;
1662                                const register double tmp2_1 = A_01_1*w4;
1663                                const register double tmp44_1 = tmp7_0*w7;
1664                                const register double tmp26_1 = tmp10_0*w4;
1665                                const register double tmp52_1 = tmp18_0*w8;
1666                                const register double tmp48_1 = A_10_1*w7;
1667                                const register double tmp46_1 = A_01_3*w8;
1668                                const register double tmp50_1 = A_01_0*w2;
1669                                const register double tmp8_1 = tmp4_0*w5;
1670                                const register double tmp56_1 = tmp19_0*w8;
1671                                const register double tmp9_1 = tmp2_0*w10;
1672                                const register double tmp19_1 = A_10_3*w2;
1673                                const register double tmp47_1 = A_10_2*w4;
1674                                const register double tmp16_1 = tmp3_0*w0;
1675                                const register double tmp18_1 = tmp1_0*w6;
1676                                const register double tmp31_1 = tmp11_0*w12;
1677                                const register double tmp55_1 = tmp15_0*w2;
1678                                const register double tmp39_1 = A_10_2*w7;
1679                                const register double tmp11_1 = tmp6_0*w7;
1680                                const register double tmp40_1 = tmp11_0*w17;
1681                                const register double tmp34_1 = tmp15_0*w8;
1682                                const register double tmp33_1 = tmp14_0*w5;
1683                                const register double tmp24_1 = tmp11_0*w13;
1684                                const register double tmp3_1 = tmp1_0*w0;
1685                                const register double tmp5_1 = tmp2_0*w3;
1686                                const register double tmp43_1 = tmp17_0*w5;
1687                                const register double tmp15_1 = A_01_2*w4;
1688                                const register double tmp53_1 = tmp19_0*w2;
1689                                const register double tmp27_1 = tmp3_0*w11;
1690                                const register double tmp32_1 = tmp13_0*w2;
1691                                const register double tmp10_1 = tmp5_0*w9;
1692                                const register double tmp37_1 = A_10_1*w4;
1693                                const register double tmp38_1 = tmp5_0*w15;
1694                                const register double tmp17_1 = A_01_1*w7;
1695                                const register double tmp12_1 = tmp7_0*w4;
1696                                const register double tmp22_1 = tmp10_0*w7;
1697                                const register double tmp57_1 = tmp18_0*w2;
1698                                const register double tmp28_1 = tmp9_0*w7;
1699                                const register double tmp29_1 = tmp1_0*w14;
1700                                const register double tmp51_1 = tmp11_0*w16;
1701                                const register double tmp42_1 = tmp12_0*w16;
1702                                const register double tmp49_1 = tmp12_0*w17;
1703                                const register double tmp21_1 = tmp1_0*w11;
1704                                const register double tmp1_1 = tmp0_0*w1;
1705                                const register double tmp45_1 = tmp6_0*w4;
1706                                const register double tmp7_1 = A_10_0*w2;
1707                                const register double tmp6_1 = tmp3_0*w6;
1708                                const register double tmp13_1 = tmp8_0*w1;
1709                                const register double tmp36_1 = tmp16_0*w1;
1710                                const register double tmp41_1 = A_01_3*w2;
1711                                const register double tmp30_1 = tmp12_0*w13;
1712                                const register double tmp4_1 = A_01_2*w7;
1713                                const register double tmp0_1 = A_10_3*w8;
1714                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
1715                                EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
1716                                EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
1717                                EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1718                                EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1719                                EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
1720                                EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
1721                                EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
1722                                EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1723                                EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
1724                                EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
1725                                EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
1726                                EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
1727                                EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
1728                                EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
1729                                EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
1730                            } else { /* constant data */
1731                                const register double A_00 = A_p[INDEX2(0,0,2)];
1732                                const register double A_01 = A_p[INDEX2(0,1,2)];
1733                                const register double A_10 = A_p[INDEX2(1,0,2)];
1734                                const register double A_11 = A_p[INDEX2(1,1,2)];
1735                                const register double tmp0_0 = A_01 + A_10;
1736                                const register double tmp0_1 = A_00*w18;
1737                                const register double tmp10_1 = A_01*w20;
1738                                const register double tmp12_1 = A_00*w26;
1739                                const register double tmp4_1 = A_00*w22;
1740                                const register double tmp8_1 = A_00*w24;
1741                                const register double tmp13_1 = A_10*w19;
1742                                const register double tmp9_1 = tmp0_0*w20;
1743                                const register double tmp3_1 = A_11*w21;
1744                                const register double tmp11_1 = A_11*w27;
1745                                const register double tmp1_1 = A_01*w19;
1746                                const register double tmp6_1 = A_11*w23;
1747                                const register double tmp7_1 = A_11*w25;
1748                                const register double tmp2_1 = A_10*w20;
1749                                const register double tmp5_1 = tmp0_0*w19;
1750                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1751                                EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1752                                EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1753                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1754                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1755                                EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1756                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1757                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1758                                EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1759                                EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1760                                EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1761                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1762                                EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1763                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1764                                EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1765                                EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1766                            }
1767                        }
1768                        ///////////////
1769                        // process B //
1770                        ///////////////
1771                        if (!B.isEmpty()) {
1772                            add_EM_S=true;
1773                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
1774                            if (B.actsExpanded()) {
1775                                const register double B_0_0 = B_p[INDEX2(0,0,2)];
1776                                const register double B_1_0 = B_p[INDEX2(1,0,2)];
1777                                const register double B_0_1 = B_p[INDEX2(0,1,2)];
1778                                const register double B_1_1 = B_p[INDEX2(1,1,2)];
1779                                const register double B_0_2 = B_p[INDEX2(0,2,2)];
1780                                const register double B_1_2 = B_p[INDEX2(1,2,2)];
1781                                const register double B_0_3 = B_p[INDEX2(0,3,2)];
1782                                const register double B_1_3 = B_p[INDEX2(1,3,2)];
1783                                const register double tmp3_0 = B_0_0 + B_0_2;
1784                                const register double tmp1_0 = B_1_2 + B_1_3;
1785                                const register double tmp2_0 = B_0_1 + B_0_3;
1786                                const register double tmp0_0 = B_1_0 + B_1_1;
1787                                const register double tmp63_1 = B_1_1*w42;
1788                                const register double tmp79_1 = B_1_1*w40;
1789                                const register double tmp37_1 = tmp3_0*w35;
1790                                const register double tmp8_1 = tmp0_0*w32;
1791                                const register double tmp71_1 = B_0_1*w34;
1792                                const register double tmp19_1 = B_0_3*w31;
1793                                const register double tmp15_1 = B_0_3*w34;
1794                                const register double tmp9_1 = tmp3_0*w34;
1795                                const register double tmp35_1 = B_1_0*w36;
1796                                const register double tmp66_1 = B_0_3*w28;
1797                                const register double tmp28_1 = B_1_0*w42;
1798                                const register double tmp22_1 = B_1_0*w40;
1799                                const register double tmp16_1 = B_1_2*w29;
1800                                const register double tmp6_1 = tmp2_0*w35;
1801                                const register double tmp55_1 = B_1_3*w40;
1802                                const register double tmp50_1 = B_1_3*w42;
1803                                const register double tmp7_1 = tmp1_0*w29;
1804                                const register double tmp1_1 = tmp1_0*w32;
1805                                const register double tmp57_1 = B_0_3*w30;
1806                                const register double tmp18_1 = B_1_1*w32;
1807                                const register double tmp53_1 = B_1_0*w41;
1808                                const register double tmp61_1 = B_1_3*w36;
1809                                const register double tmp27_1 = B_0_3*w38;
1810                                const register double tmp64_1 = B_0_2*w30;
1811                                const register double tmp76_1 = B_0_1*w38;
1812                                const register double tmp39_1 = tmp2_0*w34;
1813                                const register double tmp62_1 = B_0_1*w31;
1814                                const register double tmp56_1 = B_0_0*w31;
1815                                const register double tmp49_1 = B_1_1*w36;
1816                                const register double tmp2_1 = B_0_2*w31;
1817                                const register double tmp23_1 = B_0_2*w33;
1818                                const register double tmp38_1 = B_1_1*w43;
1819                                const register double tmp74_1 = B_1_2*w41;
1820                                const register double tmp43_1 = B_1_1*w41;
1821                                const register double tmp58_1 = B_0_2*w28;
1822                                const register double tmp67_1 = B_0_0*w33;
1823                                const register double tmp33_1 = tmp0_0*w39;
1824                                const register double tmp4_1 = B_0_0*w28;
1825                                const register double tmp20_1 = B_0_0*w30;
1826                                const register double tmp13_1 = B_0_2*w38;
1827                                const register double tmp65_1 = B_1_2*w43;
1828                                const register double tmp0_1 = tmp0_0*w29;
1829                                const register double tmp41_1 = tmp3_0*w33;
1830                                const register double tmp73_1 = B_0_2*w37;
1831                                const register double tmp69_1 = B_0_0*w38;
1832                                const register double tmp48_1 = B_1_2*w39;
1833                                const register double tmp59_1 = B_0_1*w33;
1834                                const register double tmp17_1 = B_1_3*w41;
1835                                const register double tmp5_1 = B_0_3*w33;
1836                                const register double tmp3_1 = B_0_1*w30;
1837                                const register double tmp21_1 = B_0_1*w28;
1838                                const register double tmp42_1 = B_1_0*w29;
1839                                const register double tmp54_1 = B_1_2*w32;
1840                                const register double tmp60_1 = B_1_0*w39;
1841                                const register double tmp32_1 = tmp1_0*w36;
1842                                const register double tmp10_1 = B_0_1*w37;
1843                                const register double tmp14_1 = B_0_0*w35;
1844                                const register double tmp29_1 = B_0_1*w35;
1845                                const register double tmp26_1 = B_1_2*w36;
1846                                const register double tmp30_1 = B_1_3*w43;
1847                                const register double tmp70_1 = B_0_2*w35;
1848                                const register double tmp34_1 = B_1_3*w39;
1849                                const register double tmp51_1 = B_1_0*w43;
1850                                const register double tmp31_1 = B_0_2*w34;
1851                                const register double tmp45_1 = tmp3_0*w28;
1852                                const register double tmp11_1 = tmp1_0*w39;
1853                                const register double tmp52_1 = B_1_1*w29;
1854                                const register double tmp44_1 = B_1_3*w32;
1855                                const register double tmp25_1 = B_1_1*w39;
1856                                const register double tmp47_1 = tmp2_0*w33;
1857                                const register double tmp72_1 = B_1_3*w29;
1858                                const register double tmp40_1 = tmp2_0*w28;
1859                                const register double tmp46_1 = B_1_2*w40;
1860                                const register double tmp36_1 = B_1_2*w42;
1861                                const register double tmp24_1 = B_0_0*w37;
1862                                const register double tmp77_1 = B_0_3*w35;
1863                                const register double tmp68_1 = B_0_3*w37;
1864                                const register double tmp78_1 = B_0_0*w34;
1865                                const register double tmp12_1 = tmp0_0*w36;
1866                                const register double tmp75_1 = B_1_0*w32;
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 + tmp6_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)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1877                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1878                                EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_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 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1881                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_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 B_0 = B_p[0];
1885                                const register double B_1 = B_p[1];
1886                                const register double tmp6_1 = B_1*w50;
1887                                const register double tmp1_1 = B_1*w45;
1888                                const register double tmp5_1 = B_1*w49;
1889                                const register double tmp4_1 = B_1*w48;
1890                                const register double tmp0_1 = B_0*w44;
1891                                const register double tmp2_1 = B_0*w46;
1892                                const register double tmp7_1 = B_0*w51;
1893                                const register double tmp3_1 = B_0*w47;
1894                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1895                                EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
1896                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1897                                EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1898                                EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
1899                                EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
1900                                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
1901                                EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
1902                                EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
1903                                EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
1904                                EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
1905                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
1906                                EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
1907                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1908                                EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
1909                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
1910                            }
1911                        }
1912                        ///////////////
1913                        // process C //
1914                        ///////////////
1915                        if (!C.isEmpty()) {
1916                            add_EM_S=true;
1917                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
1918                            if (C.actsExpanded()) {
1919                                const register double C_0_0 = C_p[INDEX2(0,0,2)];
1920                                const register double C_1_0 = C_p[INDEX2(1,0,2)];
1921                                const register double C_0_1 = C_p[INDEX2(0,1,2)];
1922                                const register double C_1_1 = C_p[INDEX2(1,1,2)];
1923                                const register double C_0_2 = C_p[INDEX2(0,2,2)];
1924                                const register double C_1_2 = C_p[INDEX2(1,2,2)];
1925                                const register double C_0_3 = C_p[INDEX2(0,3,2)];
1926                                const register double C_1_3 = C_p[INDEX2(1,3,2)];
1927                                const register double tmp2_0 = C_0_1 + C_0_3;
1928                                const register double tmp1_0 = C_1_2 + C_1_3;
1929                                const register double tmp3_0 = C_0_0 + C_0_2;
1930                                const register double tmp0_0 = C_1_0 + C_1_1;
1931                                const register double tmp64_1 = C_0_2*w30;
1932                                const register double tmp14_1 = C_0_2*w28;
1933                                const register double tmp19_1 = C_0_3*w31;
1934                                const register double tmp22_1 = C_1_0*w40;
1935                                const register double tmp37_1 = tmp3_0*w35;
1936                                const register double tmp29_1 = C_0_1*w35;
1937                                const register double tmp73_1 = C_0_2*w37;
1938                                const register double tmp74_1 = C_1_2*w41;
1939                                const register double tmp52_1 = C_1_3*w39;
1940                                const register double tmp25_1 = C_1_1*w39;
1941                                const register double tmp62_1 = C_0_1*w31;
1942                                const register double tmp79_1 = C_1_1*w40;
1943                                const register double tmp43_1 = C_1_1*w36;
1944                                const register double tmp27_1 = C_0_3*w38;
1945                                const register double tmp28_1 = C_1_0*w42;
1946                                const register double tmp63_1 = C_1_1*w42;
1947                                const register double tmp59_1 = C_0_3*w34;
1948                                const register double tmp72_1 = C_1_3*w29;
1949                                const register double tmp40_1 = tmp2_0*w35;
1950                                const register double tmp13_1 = C_0_3*w30;
1951                                const register double tmp51_1 = C_1_2*w40;
1952                                const register double tmp54_1 = C_1_2*w42;
1953                                const register double tmp12_1 = C_0_0*w31;
1954                                const register double tmp2_1 = tmp1_0*w32;
1955                                const register double tmp68_1 = C_0_2*w31;
1956                                const register double tmp75_1 = C_1_0*w32;
1957                                const register double tmp49_1 = C_1_1*w41;
1958                                const register double tmp4_1 = C_0_2*w35;
1959                                const register double tmp66_1 = C_0_3*w28;
1960                                const register double tmp56_1 = C_0_1*w37;
1961                                const register double tmp5_1 = C_0_1*w34;
1962                                const register double tmp38_1 = tmp2_0*w34;
1963                                const register double tmp76_1 = C_0_1*w38;
1964                                const register double tmp21_1 = C_0_1*w28;
1965                                const register double tmp69_1 = C_0_1*w30;
1966                                const register double tmp53_1 = C_1_0*w36;
1967                                const register double tmp42_1 = C_1_2*w39;
1968                                const register double tmp32_1 = tmp1_0*w29;
1969                                const register double tmp45_1 = C_1_0*w43;
1970                                const register double tmp33_1 = tmp0_0*w32;
1971                                const register double tmp35_1 = C_1_0*w41;
1972                                const register double tmp26_1 = C_1_2*w36;
1973                                const register double tmp67_1 = C_0_0*w33;
1974                                const register double tmp31_1 = C_0_2*w34;
1975                                const register double tmp20_1 = C_0_0*w30;
1976                                const register double tmp70_1 = C_0_0*w28;
1977                                const register double tmp8_1 = tmp0_0*w39;
1978                                const register double tmp30_1 = C_1_3*w43;
1979                                const register double tmp0_1 = tmp0_0*w29;
1980                                const register double tmp17_1 = C_1_3*w41;
1981                                const register double tmp58_1 = C_0_0*w35;
1982                                const register double tmp9_1 = tmp3_0*w33;
1983                                const register double tmp61_1 = C_1_3*w36;
1984                                const register double tmp41_1 = tmp3_0*w34;
1985                                const register double tmp50_1 = C_1_3*w32;
1986                                const register double tmp18_1 = C_1_1*w32;
1987                                const register double tmp6_1 = tmp1_0*w36;
1988                                const register double tmp3_1 = C_0_0*w38;
1989                                const register double tmp34_1 = C_1_1*w29;
1990                                const register double tmp77_1 = C_0_3*w35;
1991                                const register double tmp65_1 = C_1_2*w43;
1992                                const register double tmp71_1 = C_0_3*w33;
1993                                const register double tmp55_1 = C_1_1*w43;
1994                                const register double tmp46_1 = tmp3_0*w28;
1995                                const register double tmp24_1 = C_0_0*w37;
1996                                const register double tmp10_1 = tmp1_0*w39;
1997                                const register double tmp48_1 = C_1_0*w29;
1998                                const register double tmp15_1 = C_0_1*w33;
1999                                const register double tmp36_1 = C_1_2*w32;
2000                                const register double tmp60_1 = C_1_0*w39;
2001                                const register double tmp47_1 = tmp2_0*w33;
2002                                const register double tmp16_1 = C_1_2*w29;
2003                                const register double tmp1_1 = C_0_3*w37;
2004                                const register double tmp7_1 = tmp2_0*w28;
2005                                const register double tmp39_1 = C_1_3*w40;
2006                                const register double tmp44_1 = C_1_3*w42;
2007                                const register double tmp57_1 = C_0_2*w38;
2008                                const register double tmp78_1 = C_0_0*w34;
2009                                const register double tmp11_1 = tmp0_0*w36;
2010                                const register double tmp23_1 = C_0_2*w33;
2011                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
2012                                EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
2013                                EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2014                                EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
2015                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2016                                EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
2017                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
2018                                EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
2019                                EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
2020                                EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2021                                EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
2022                                EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
2023                                EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
2024                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
2025                                EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
2026                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
2027                            } else { /* constant data */
2028                                const register double C_0 = C_p[0];
2029                                const register double C_1 = C_p[1];
2030                                const register double tmp1_1 = C_1*w45;
2031                                const register double tmp3_1 = C_0*w51;
2032                                const register double tmp4_1 = C_0*w44;
2033                                const register double tmp7_1 = C_0*w46;
2034                                const register double tmp5_1 = C_1*w49;
2035                                const register double tmp2_1 = C_1*w48;
2036                                const register double tmp0_1 = C_0*w47;
2037                                const register double tmp6_1 = C_1*w50;
2038                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2039                                EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
2040                                EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
2041                                EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
2042                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
2043                                EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
2044                                EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
2045                                EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
2046                                EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
2047                                EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
2048                                EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
2049                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2050                                EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
2051                                EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
2052                                EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
2053                                EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
2054                            }
2055                        }
2056                        ///////////////
2057                        // process D //
2058                        ///////////////
2059                        if (!D.isEmpty()) {
2060                            add_EM_S=true;
2061                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2062                            if (D.actsExpanded()) {
2063                                const register double D_0 = D_p[0];
2064                                const register double D_1 = D_p[1];
2065                                const register double D_2 = D_p[2];
2066                                const register double D_3 = D_p[3];
2067                                const register double tmp4_0 = D_1 + D_3;
2068                                const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;
2069                                const register double tmp5_0 = D_0 + D_2;
2070                                const register double tmp0_0 = D_0 + D_1;
2071                                const register double tmp6_0 = D_0 + D_3;
2072                                const register double tmp1_0 = D_2 + D_3;
2073                                const register double tmp3_0 = D_1 + D_2;
2074                                const register double tmp16_1 = D_1*w56;
2075                                const register double tmp14_1 = tmp6_0*w54;
2076                                const register double tmp8_1 = D_3*w55;
2077                                const register double tmp2_1 = tmp2_0*w54;
2078                                const register double tmp12_1 = tmp5_0*w52;
2079                                const register double tmp4_1 = tmp0_0*w53;
2080                                const register double tmp3_1 = tmp1_0*w52;
2081                                const register double tmp13_1 = tmp4_0*w53;
2082                                const register double tmp10_1 = tmp4_0*w52;
2083                                const register double tmp1_1 = tmp1_0*w53;
2084                                const register double tmp7_1 = D_3*w56;
2085                                const register double tmp0_1 = tmp0_0*w52;
2086                                const register double tmp11_1 = tmp5_0*w53;
2087                                const register double tmp9_1 = D_0*w56;
2088                                const register double tmp5_1 = tmp3_0*w54;
2089                                const register double tmp18_1 = D_2*w56;
2090                                const register double tmp17_1 = D_1*w55;
2091                                const register double tmp6_1 = D_0*w55;
2092                                const register double tmp15_1 = D_2*w55;
2093                                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2094                                EM_S[INDEX2(1,2,4)]+=tmp2_1;
2095                                EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
2096                                EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2097                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2098                                EM_S[INDEX2(3,0,4)]+=tmp2_1;
2099                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
2100                                EM_S[INDEX2(2,1,4)]+=tmp2_1;
2101                                EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
2102                                EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
2103                                EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
2104                                EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2105                                EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
2106                                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2107                                EM_S[INDEX2(0,3,4)]+=tmp2_1;
2108                                EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
2109                            } else { /* constant data */
2110                                const register double D_0 = D_p[0];
2111                                const register double tmp2_1 = D_0*w59;
2112                                const register double tmp1_1 = D_0*w58;
2113                                const register double tmp0_1 = D_0*w57;
2114                                EM_S[INDEX2(0,1,4)]+=tmp0_1;
2115                                EM_S[INDEX2(1,2,4)]+=tmp1_1;
2116                                EM_S[INDEX2(3,2,4)]+=tmp0_1;
2117                                EM_S[INDEX2(0,0,4)]+=tmp2_1;
2118                                EM_S[INDEX2(3,3,4)]+=tmp2_1;
2119                                EM_S[INDEX2(3,0,4)]+=tmp1_1;
2120                                EM_S[INDEX2(3,1,4)]+=tmp0_1;
2121                                EM_S[INDEX2(2,1,4)]+=tmp1_1;
2122                                EM_S[INDEX2(0,2,4)]+=tmp0_1;
2123                                EM_S[INDEX2(2,0,4)]+=tmp0_1;
2124                                EM_S[INDEX2(1,3,4)]+=tmp0_1;
2125                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
2126                                EM_S[INDEX2(2,2,4)]+=tmp2_1;
2127                                EM_S[INDEX2(1,0,4)]+=tmp0_1;
2128                                EM_S[INDEX2(0,3,4)]+=tmp1_1;
2129                                EM_S[INDEX2(1,1,4)]+=tmp2_1;
2130                            }
2131                        }
2132                        ///////////////
2133                        // process X //
2134                        ///////////////
2135                        if (!X.isEmpty()) {
2136                            add_EM_F=true;
2137                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2138                            if (X.actsExpanded()) {
2139                                const register double X_0_0 = X_p[INDEX2(0,0,2)];
2140                                const register double X_1_0 = X_p[INDEX2(1,0,2)];
2141                                const register double X_0_1 = X_p[INDEX2(0,1,2)];
2142                                const register double X_1_1 = X_p[INDEX2(1,1,2)];
2143                                const register double X_0_2 = X_p[INDEX2(0,2,2)];
2144                                const register double X_1_2 = X_p[INDEX2(1,2,2)];
2145                                const register double X_0_3 = X_p[INDEX2(0,3,2)];
2146                                const register double X_1_3 = X_p[INDEX2(1,3,2)];
2147                                const register double tmp1_0 = X_1_1 + X_1_3;
2148                                const register double tmp3_0 = X_0_0 + X_0_1;
2149                                const register double tmp2_0 = X_1_0 + X_1_2;
2150                                const register double tmp0_0 = X_0_2 + X_0_3;
2151                                const register double tmp8_1 = tmp2_0*w66;
2152                                const register double tmp5_1 = tmp3_0*w64;
2153                                const register double tmp14_1 = tmp0_0*w64;
2154                                const register double tmp3_1 = tmp3_0*w60;
2155                                const register double tmp9_1 = tmp3_0*w63;
2156                                const register double tmp13_1 = tmp3_0*w65;
2157                                const register double tmp12_1 = tmp1_0*w66;
2158                                const register double tmp10_1 = tmp0_0*w60;
2159                                const register double tmp2_1 = tmp2_0*w61;
2160                                const register double tmp6_1 = tmp2_0*w62;
2161                                const register double tmp4_1 = tmp0_0*w65;
2162                                const register double tmp11_1 = tmp1_0*w67;
2163                                const register double tmp1_1 = tmp1_0*w62;
2164                                const register double tmp7_1 = tmp1_0*w61;
2165                                const register double tmp0_1 = tmp0_0*w63;
2166                                const register double tmp15_1 = tmp2_0*w67;
2167                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2168                                EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2169                                EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2170                                EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2171                            } else { /* constant data */
2172                                const register double X_0 = X_p[0];
2173                                const register double X_1 = X_p[1];
2174                                const register double tmp3_1 = X_1*w71;
2175                                const register double tmp2_1 = X_0*w70;
2176                                const register double tmp1_1 = X_0*w68;
2177                                const register double tmp0_1 = X_1*w69;
2178                                EM_F[0]+=tmp0_1 + tmp1_1;
2179                                EM_F[1]+=tmp0_1 + tmp2_1;
2180                                EM_F[2]+=tmp1_1 + tmp3_1;
2181                                EM_F[3]+=tmp2_1 + tmp3_1;
2182                            }
2183                        }
2184                        ///////////////
2185                        // process Y //
2186                        ///////////////
2187                        if (!Y.isEmpty()) {
2188                            add_EM_F=true;
2189                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2190                            if (Y.actsExpanded()) {
2191                                const register double Y_0 = Y_p[0];
2192                                const register double Y_1 = Y_p[1];
2193                                const register double Y_2 = Y_p[2];
2194                                const register double Y_3 = Y_p[3];
2195                                const register double tmp0_0 = Y_1 + Y_2;
2196                                const register double tmp1_0 = Y_0 + Y_3;
2197                                const register double tmp9_1 = Y_0*w74;
2198                                const register double tmp4_1 = tmp1_0*w73;
2199                                const register double tmp5_1 = Y_2*w74;
2200                                const register double tmp7_1 = Y_1*w74;
2201                                const register double tmp6_1 = Y_2*w72;
2202                                const register double tmp2_1 = Y_3*w74;
2203                                const register double tmp8_1 = Y_3*w72;
2204                                const register double tmp3_1 = Y_1*w72;
2205                                const register double tmp0_1 = Y_0*w72;
2206                                const register double tmp1_1 = tmp0_0*w73;
2207                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2208                                EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2209                                EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2210                                EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2211                            } else { /* constant data */
2212                                const register double Y_0 = Y_p[0];
2213                                const register double tmp0_1 = Y_0*w75;
2214                                EM_F[0]+=tmp0_1;
2215                                EM_F[1]+=tmp0_1;
2216                                EM_F[2]+=tmp0_1;
2217                                EM_F[3]+=tmp0_1;
2218                            }
2219                        }
2220                        /* GENERATOR SNIP_PDE_SINGLE BOTTOM */
2221    
2222                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2223                        IndexVector rowIndex;
2224                        const index_t firstNode=k1*m_N0+k0;
2225                        rowIndex.push_back(firstNode);
2226                        rowIndex.push_back(firstNode+1);
2227                        rowIndex.push_back(firstNode+m_N0);
2228                        rowIndex.push_back(firstNode+1+m_N0);
2229                        if (add_EM_S) {
2230                            addToSystemMatrix(mat, 4, rowIndex, 1, 4, rowIndex, 1, EM_S);
2231                        }
2232                        if (add_EM_F) {
2233                            double *F_p=rhs.getSampleDataRW(0);
2234                            for (index_t i=0; i<4; i++) {
2235                                if (rowIndex[i]<getNumNodes())
2236                                    F_p[rowIndex[i]]+=EM_F[i];
2237                            }
2238                        }
2239                    } // end k0 loop
2240                } // end k1 loop
2241            } // end of coloring
2242        } // end of parallel region
2243    }
2244    
2245    void Rectangle::addToSystemMatrix(Paso_SystemMatrix* in, dim_t NN_Eq,
2246           const IndexVector& nodes_Eq, dim_t num_Eq, dim_t NN_Sol,
2247           const IndexVector& nodes_Sol, dim_t num_Sol, const vector<double>& array) const
2248    {
2249        const dim_t row_block_size = in->row_block_size;
2250        const dim_t col_block_size = in->col_block_size;
2251        const dim_t block_size = in->block_size;
2252        const dim_t num_subblocks_Eq = num_Eq / row_block_size;
2253        const dim_t num_subblocks_Sol = num_Sol / col_block_size;
2254        const dim_t numMyCols = in->pattern->mainPattern->numInput;
2255        const dim_t numMyRows = in->pattern->mainPattern->numOutput;
2256    
2257        const index_t* mainBlock_ptr = in->mainBlock->pattern->ptr;
2258        const index_t* mainBlock_index = in->mainBlock->pattern->index;
2259        double* mainBlock_val = in->mainBlock->val;
2260        const index_t* col_coupleBlock_ptr = in->col_coupleBlock->pattern->ptr;
2261        const index_t* col_coupleBlock_index = in->col_coupleBlock->pattern->index;
2262        double* col_coupleBlock_val = in->col_coupleBlock->val;
2263        const index_t* row_coupleBlock_ptr = in->row_coupleBlock->pattern->ptr;
2264        const index_t* row_coupleBlock_index = in->row_coupleBlock->pattern->index;
2265        double* row_coupleBlock_val = in->row_coupleBlock->val;
2266        for (dim_t k_Eq = 0; k_Eq < NN_Eq; ++k_Eq) {
2267            // down columns of array
2268            const dim_t j_Eq = nodes_Eq[k_Eq];
2269            for (dim_t l_row = 0; l_row < num_subblocks_Eq; ++l_row) {
2270                const dim_t i_row = j_Eq * num_subblocks_Eq + l_row;
2271                // only look at the matrix rows stored on this processor
2272                if (i_row < numMyRows) {
2273                    for (dim_t k_Sol = 0; k_Sol < NN_Sol; ++k_Sol) {
2274                        // across rows of array
2275                        const dim_t j_Sol = nodes_Sol[k_Sol];
2276                        for (dim_t l_col = 0; l_col < num_subblocks_Sol; ++l_col) {
2277                            // only look at the matrix rows stored on this processor
2278                            const dim_t i_col = j_Sol * num_subblocks_Sol + l_col;
2279                            if (i_col < numMyCols) {
2280                                for (dim_t k = mainBlock_ptr[i_row];
2281                                     k < mainBlock_ptr[i_row + 1]; ++k) {
2282                                    if (mainBlock_index[k] == i_col) {
2283                                        // entry array(k_Sol, j_Eq) is a block
2284                                        // (row_block_size x col_block_size)
2285                                        for (dim_t ic = 0; ic < col_block_size; ++ic) {
2286                                            const dim_t i_Sol = ic + col_block_size * l_col;
2287                                            for (dim_t ir = 0; ir < row_block_size; ++ir) {
2288                                                const dim_t i_Eq = ir + row_block_size * l_row;
2289                                                mainBlock_val[k*block_size + ir + row_block_size*ic] +=
2290                                                    array[INDEX4
2291                                                          (i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, NN_Eq)];
2292                                            }
2293                                        }
2294                                        break;
2295                                    }
2296                                }
2297                            } else {
2298                                for (dim_t k = col_coupleBlock_ptr[i_row];
2299                                     k < col_coupleBlock_ptr[i_row + 1]; ++k) {
2300                                    if (col_coupleBlock_index[k] == i_col - numMyCols) {
2301                                        // entry array(k_Sol, j_Eq) is a block
2302                                        // (row_block_size x col_block_size)
2303                                        for (dim_t ic = 0; ic < col_block_size; ++ic) {
2304                                            const dim_t i_Sol = ic + col_block_size * l_col;
2305                                            for (dim_t ir = 0; ir < row_block_size; ++ir) {
2306                                                const dim_t i_Eq = ir + row_block_size * l_row;
2307                                                col_coupleBlock_val[k * block_size + ir + row_block_size * ic] +=
2308                                                    array[INDEX4
2309                                                          (i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, NN_Eq)];
2310                                            }
2311                                        }
2312                                        break;
2313                                    }
2314                                }
2315                            }
2316                        }
2317                    }
2318                } else {
2319                    for (dim_t k_Sol = 0; k_Sol < NN_Sol; ++k_Sol) {
2320                        // across rows of array
2321                        const dim_t j_Sol = nodes_Sol[k_Sol];
2322                        for (dim_t l_col = 0; l_col < num_subblocks_Sol; ++l_col) {
2323                            const dim_t i_col = j_Sol * num_subblocks_Sol + l_col;
2324                            if (i_col < numMyCols) {
2325                                for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
2326                                     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
2327                                {
2328                                    if (row_coupleBlock_index[k] == i_col) {
2329                                        // entry array(k_Sol, j_Eq) is a block
2330                                        // (row_block_size x col_block_size)
2331                                        for (dim_t ic = 0; ic < col_block_size; ++ic) {
2332                                            const dim_t i_Sol = ic + col_block_size * l_col;
2333                                            for (dim_t ir = 0; ir < row_block_size; ++ir) {
2334                                                const dim_t i_Eq = ir + row_block_size * l_row;
2335                                                row_coupleBlock_val[k * block_size + ir + row_block_size * ic] +=
2336                                                    array[INDEX4
2337                                                          (i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, NN_Eq)];
2338                                            }
2339                                        }
2340                                        break;
2341                                    }
2342                                }
2343                            }
2344                        }
2345                    }
2346                }
2347            }
2348        }
2349  }  }
2350    
2351  } // end of namespace ripley  } // end of namespace ripley

Legend:
Removed from v.3697  
changed lines
  Added in v.3748

  ViewVC Help
Powered by ViewVC 1.1.26