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

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

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

revision 3762 by caltinay, Tue Jan 10 00:01:45 2012 UTC revision 3764 by caltinay, Tue Jan 10 06:31:15 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2011 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 49  Rectangle::Rectangle(int n0, int n1, dou Line 49  Rectangle::Rectangle(int n0, int n1, dou
49      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
50          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
51    
52      // local number of elements (including overlap)      // local number of elements (with and without overlap)
53      m_NE0 = (m_NX>1 ? (n0+1)/m_NX : n0);      m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
54      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)      if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
55          m_NE0++;          m_NE0++;
56      m_NE1 = (m_NY>1 ? (n1+1)/m_NY : n1);      else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)
57            m_ownNE0--;
58    
59        m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
60      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)      if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)
61          m_NE1++;          m_NE1++;
62        else if (m_NY>1 && m_mpiInfo->rank/m_NX==m_NY-1)
63            m_ownNE1--;
64    
65      // local number of nodes      // local number of nodes
66      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
# Line 281  bool Rectangle::ownSample(int fsType, in Line 286  bool Rectangle::ownSample(int fsType, in
286          case FaceElements:          case FaceElements:
287          case ReducedFaceElements:          case ReducedFaceElements:
288              {              {
289                  // check ownership of face element's first node                  // determine which face the sample belongs to before
290                    // checking ownership of face element's first node
291                  const IndexVector faces = getNumFacesPerBoundary();                  const IndexVector faces = getNumFacesPerBoundary();
292                  dim_t n=0;                  dim_t n=0;
293                  for (size_t i=0; i<faces.size(); i++) {                  for (size_t i=0; i<faces.size(); i++) {
# Line 311  bool Rectangle::ownSample(int fsType, in Line 317  bool Rectangle::ownSample(int fsType, in
317      throw RipleyException(msg.str());      throw RipleyException(msg.str());
318  }  }
319    
320  void Rectangle::setToGradient(escript::Data& out, const escript::Data& cIn) const  void Rectangle::setToNormal(escript::Data& out) const
321    {
322        if (out.getFunctionSpace().getTypeCode() == FaceElements) {
323            out.requireWrite();
324    #pragma omp parallel
325            {
326                if (m_faceOffset[0] > -1) {
327    #pragma omp for nowait
328                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
329                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
330                        // set vector at two quadrature points
331                        *o++ = -1.;
332                        *o++ = 0.;
333                        *o++ = -1.;
334                        *o = 0.;
335                    }
336                }
337    
338                if (m_faceOffset[1] > -1) {
339    #pragma omp for nowait
340                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
341                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
342                        // set vector at two quadrature points
343                        *o++ = 1.;
344                        *o++ = 0.;
345                        *o++ = 1.;
346                        *o = 0.;
347                    }
348                }
349    
350                if (m_faceOffset[2] > -1) {
351    #pragma omp for nowait
352                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
353                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
354                        // set vector at two quadrature points
355                        *o++ = 0.;
356                        *o++ = -1.;
357                        *o++ = 0.;
358                        *o = -1.;
359                    }
360                }
361    
362                if (m_faceOffset[3] > -1) {
363    #pragma omp for nowait
364                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
365                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
366                        // set vector at two quadrature points
367                        *o++ = 0.;
368                        *o++ = 1.;
369                        *o++ = 0.;
370                        *o = 1.;
371                    }
372                }
373            } // end of parallel section
374        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
375            out.requireWrite();
376    #pragma omp parallel
377            {
378                if (m_faceOffset[0] > -1) {
379    #pragma omp for nowait
380                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
381                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
382                        *o++ = -1.;
383                        *o = 0.;
384                    }
385                }
386    
387                if (m_faceOffset[1] > -1) {
388    #pragma omp for nowait
389                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
390                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
391                        *o++ = 1.;
392                        *o = 0.;
393                    }
394                }
395    
396                if (m_faceOffset[2] > -1) {
397    #pragma omp for nowait
398                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
399                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
400                        *o++ = 0.;
401                        *o = -1.;
402                    }
403                }
404    
405                if (m_faceOffset[3] > -1) {
406    #pragma omp for nowait
407                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
408                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
409                        *o++ = 0.;
410                        *o = 1.;
411                    }
412                }
413            } // end of parallel section
414    
415        } else {
416            stringstream msg;
417            msg << "setToNormal() not implemented for "
418                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
419            throw RipleyException(msg.str());
420        }
421    }
422    
423    void Rectangle::setToSize(escript::Data& out) const
424    {
425        if (out.getFunctionSpace().getTypeCode() == Elements
426                || out.getFunctionSpace().getTypeCode() == ReducedElements) {
427            out.requireWrite();
428            const dim_t numQuad=out.getNumDataPointsPerSample();
429            const double hSize=getFirstCoordAndSpacing(0).second;
430            const double vSize=getFirstCoordAndSpacing(1).second;
431            const double size=min(hSize,vSize);
432    #pragma omp parallel for
433            for (index_t k = 0; k < getNumElements(); ++k) {
434                double* o = out.getSampleDataRW(k);
435                fill(o, o+numQuad, size);
436            }
437        } else if (out.getFunctionSpace().getTypeCode() == FaceElements
438                || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
439            out.requireWrite();
440            const dim_t numQuad=out.getNumDataPointsPerSample();
441            const double hSize=getFirstCoordAndSpacing(0).second;
442            const double vSize=getFirstCoordAndSpacing(1).second;
443    #pragma omp parallel
444            {
445                if (m_faceOffset[0] > -1) {
446    #pragma omp for nowait
447                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
448                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
449                        fill(o, o+numQuad, vSize);
450                    }
451                }
452    
453                if (m_faceOffset[1] > -1) {
454    #pragma omp for nowait
455                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
456                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
457                        fill(o, o+numQuad, vSize);
458                    }
459                }
460    
461                if (m_faceOffset[2] > -1) {
462    #pragma omp for nowait
463                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
464                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
465                        fill(o, o+numQuad, hSize);
466                    }
467                }
468    
469                if (m_faceOffset[3] > -1) {
470    #pragma omp for nowait
471                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
472                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
473                        fill(o, o+numQuad, hSize);
474                    }
475                }
476            } // end of parallel section
477    
478        } else {
479            stringstream msg;
480            msg << "setToSize() not implemented for "
481                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
482            throw RipleyException(msg.str());
483        }
484    }
485    
486    Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
487                                                    bool reducedColOrder) const
488    {
489        if (reducedRowOrder || reducedColOrder)
490            throw RipleyException("getPattern() not implemented for reduced order");
491    
492        return m_pattern;
493    }
494    
495    void Rectangle::Print_Mesh_Info(const bool full) const
496    {
497        RipleyDomain::Print_Mesh_Info(full);
498        if (full) {
499            cout << "     Id  Coordinates" << endl;
500            cout.precision(15);
501            cout.setf(ios::scientific, ios::floatfield);
502            pair<double,double> xdx = getFirstCoordAndSpacing(0);
503            pair<double,double> ydy = getFirstCoordAndSpacing(1);
504            for (index_t i=0; i < getNumNodes(); i++) {
505                cout << "  " << setw(5) << m_nodeId[i]
506                    << "  " << xdx.first+(i%m_N0)*xdx.second
507                    << "  " << ydy.first+(i/m_N0)*ydy.second << endl;
508            }
509        }
510    }
511    
512    IndexVector Rectangle::getNumNodesPerDim() const
513    {
514        IndexVector ret;
515        ret.push_back(m_N0);
516        ret.push_back(m_N1);
517        return ret;
518    }
519    
520    IndexVector Rectangle::getNumElementsPerDim() const
521    {
522        IndexVector ret;
523        ret.push_back(m_NE0);
524        ret.push_back(m_NE1);
525        return ret;
526    }
527    
528    IndexVector Rectangle::getNumFacesPerBoundary() const
529    {
530        IndexVector ret(4, 0);
531        //left
532        if (m_offset0==0)
533            ret[0]=m_NE1;
534        //right
535        if (m_mpiInfo->rank%m_NX==m_NX-1)
536            ret[1]=m_NE1;
537        //bottom
538        if (m_offset1==0)
539            ret[2]=m_NE0;
540        //top
541        if (m_mpiInfo->rank/m_NX==m_NY-1)
542            ret[3]=m_NE0;
543        return ret;
544    }
545    
546    pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
547    {
548        if (dim==0) {
549            return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
550        } else if (dim==1) {
551            return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
552        }
553        throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
554    }
555    
556    //protected
557    dim_t Rectangle::getNumDOF() const
558    {
559        return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;
560    }
561    
562    //protected
563    dim_t Rectangle::getNumFaceElements() const
564    {
565        const IndexVector faces = getNumFacesPerBoundary();
566        dim_t n=0;
567        for (size_t i=0; i<faces.size(); i++)
568            n+=faces[i];
569        return n;
570    }
571    
572    //protected
573    void Rectangle::assembleCoordinates(escript::Data& arg) const
574    {
575        escriptDataC x = arg.getDataC();
576        int numDim = m_numDim;
577        if (!isDataPointShapeEqual(&x, 1, &numDim))
578            throw RipleyException("setToX: Invalid Data object shape");
579        if (!numSamplesEqual(&x, 1, getNumNodes()))
580            throw RipleyException("setToX: Illegal number of samples in Data object");
581    
582        pair<double,double> xdx = getFirstCoordAndSpacing(0);
583        pair<double,double> ydy = getFirstCoordAndSpacing(1);
584        arg.requireWrite();
585    #pragma omp parallel for
586        for (dim_t i1 = 0; i1 < m_N1; i1++) {
587            for (dim_t i0 = 0; i0 < m_N0; i0++) {
588                double* point = arg.getSampleDataRW(i0+m_N0*i1);
589                point[0] = xdx.first+i0*xdx.second;
590                point[1] = ydy.first+i1*ydy.second;
591            }
592        }
593    }
594    
595    //protected
596    void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const
597  {  {
     escript::Data& in = *const_cast<escript::Data*>(&cIn);  
598      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
599      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
600      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 336  void Rectangle::setToGradient(escript::D Line 617  void Rectangle::setToGradient(escript::D
617    
618      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
619          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */  
620  #pragma omp parallel for  #pragma omp parallel for
621          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
622              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 357  void Rectangle::setToGradient(escript::D Line 637  void Rectangle::setToGradient(escript::D
637                  } /* end of component loop i */                  } /* end of component loop i */
638              } /* end of k0 loop */              } /* end of k0 loop */
639          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */  
640      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
641          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */  
642  #pragma omp parallel for  #pragma omp parallel for
643          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
644              for (index_t k0=0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 375  void Rectangle::setToGradient(escript::D Line 653  void Rectangle::setToGradient(escript::D
653                  } /* end of component loop i */                  } /* end of component loop i */
654              } /* end of k0 loop */              } /* end of k0 loop */
655          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */  
656      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
657          out.requireWrite();          out.requireWrite();
658  #pragma omp parallel  #pragma omp parallel
659          {          {
             /*** GENERATOR SNIP_GRAD_FACES TOP */  
660              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
661  #pragma omp for nowait  #pragma omp for nowait
662                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 445  void Rectangle::setToGradient(escript::D Line 721  void Rectangle::setToGradient(escript::D
721                      } /* end of component loop i */                      } /* end of component loop i */
722                  } /* end of k0 loop */                  } /* end of k0 loop */
723              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_GRAD_FACES BOTTOM */  
724          } // end of parallel section          } // end of parallel section
725      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
726          out.requireWrite();          out.requireWrite();
727  #pragma omp parallel  #pragma omp parallel
728          {          {
             /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */  
729              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
730  #pragma omp for nowait  #pragma omp for nowait
731                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 508  void Rectangle::setToGradient(escript::D Line 782  void Rectangle::setToGradient(escript::D
782                      } /* end of component loop i */                      } /* end of component loop i */
783                  } /* end of k0 loop */                  } /* end of k0 loop */
784              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */  
785          } // end of parallel section          } // end of parallel section
     } else {  
         stringstream msg;  
         msg << "setToGradient() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
786      }      }
787  }  }
788    
789  void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  //protected
790    void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
791  {  {
792      escript::Data& in = *const_cast<escript::Data*>(&arg);      const dim_t numComp = arg.getDataPointSize();
     const dim_t numComp = in.getDataPointSize();  
793      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
794      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
795        const index_t left = (m_offset0==0 ? 0 : 1);
796        const index_t bottom = (m_offset1==0 ? 0 : 1);
797      if (arg.getFunctionSpace().getTypeCode() == Elements) {      if (arg.getFunctionSpace().getTypeCode() == Elements) {
798          const double w_0 = h0*h1/4.;          const double w = h0*h1/4.;
799  #pragma omp parallel  #pragma omp parallel
800          {          {
801              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
802  #pragma omp for nowait  #pragma omp for nowait
803              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
804                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
805                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
806                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
807                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
808                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
809                          const register double f_2 = f[INDEX2(i,2,numComp)];                          const register double f2 = f[INDEX2(i,2,numComp)];
810                          const register double f_3 = f[INDEX2(i,3,numComp)];                          const register double f3 = f[INDEX2(i,3,numComp)];
811                          int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;                          int_local[i]+=(f0+f1+f2+f3)*w;
812                      }  /* end of component loop i */                      }  /* end of component loop i */
813                  } /* end of k0 loop */                  } /* end of k0 loop */
814              } /* end of k1 loop */              } /* end of k1 loop */
# Line 548  void Rectangle::setToIntegrals(vector<do Line 818  void Rectangle::setToIntegrals(vector<do
818                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
819          } // end of parallel section          } // end of parallel section
820      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
821          const double w_0 = h0*h1;          const double w = h0*h1;
822  #pragma omp parallel  #pragma omp parallel
823          {          {
824              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
825  #pragma omp for nowait  #pragma omp for nowait
826              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
827                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
828                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
829                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
830                          int_local[i]+=f[i]*w_0;                          int_local[i]+=f[i]*w;
831                      }  /* end of component loop i */                      }  /* end of component loop i */
832                  } /* end of k0 loop */                  } /* end of k0 loop */
833              } /* end of k1 loop */              } /* end of k1 loop */
# Line 567  void Rectangle::setToIntegrals(vector<do Line 837  void Rectangle::setToIntegrals(vector<do
837                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
838          } // end of parallel section          } // end of parallel section
839      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
840          const double w_0 = h0/2.;          const double w0 = h0/2.;
841          const double w_1 = h1/2.;          const double w1 = h1/2.;
842  #pragma omp parallel  #pragma omp parallel
843          {          {
844              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
845              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
846  #pragma omp for nowait  #pragma omp for nowait
847                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
848                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
849                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
850                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
851                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
852                          int_local[i]+=(f_0+f_1)*w_1;                          int_local[i]+=(f0+f1)*w1;
853                      }  /* end of component loop i */                      }  /* end of component loop i */
854                  } /* end of k1 loop */                  } /* end of k1 loop */
855              }              }
856    
857              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
858  #pragma omp for nowait  #pragma omp for nowait
859                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
860                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
861                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
862                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
863                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
864                          int_local[i]+=(f_0+f_1)*w_1;                          int_local[i]+=(f0+f1)*w1;
865                      }  /* end of component loop i */                      }  /* end of component loop i */
866                  } /* end of k1 loop */                  } /* end of k1 loop */
867              }              }
868    
869              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
870  #pragma omp for nowait  #pragma omp for nowait
871                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
872                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
873                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
874                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
875                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
876                          int_local[i]+=(f_0+f_1)*w_0;                          int_local[i]+=(f0+f1)*w0;
877                      }  /* end of component loop i */                      }  /* end of component loop i */
878                  } /* end of k0 loop */                  } /* end of k0 loop */
879              }              }
880    
881              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
882  #pragma omp for nowait  #pragma omp for nowait
883                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
884                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
885                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
886                          const register double f_0 = f[INDEX2(i,0,numComp)];                          const register double f0 = f[INDEX2(i,0,numComp)];
887                          const register double f_1 = f[INDEX2(i,1,numComp)];                          const register double f1 = f[INDEX2(i,1,numComp)];
888                          int_local[i]+=(f_0+f_1)*w_0;                          int_local[i]+=(f0+f1)*w0;
889                      }  /* end of component loop i */                      }  /* end of component loop i */
890                  } /* end of k0 loop */                  } /* end of k0 loop */
891              }              }
# Line 630  void Rectangle::setToIntegrals(vector<do Line 900  void Rectangle::setToIntegrals(vector<do
900              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
901              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
902  #pragma omp for nowait  #pragma omp for nowait
903                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
904                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
905                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
906                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
907                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 640  void Rectangle::setToIntegrals(vector<do Line 910  void Rectangle::setToIntegrals(vector<do
910    
911              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
912  #pragma omp for nowait  #pragma omp for nowait
913                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
914                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
915                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
916                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
917                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 650  void Rectangle::setToIntegrals(vector<do Line 920  void Rectangle::setToIntegrals(vector<do
920    
921              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
922  #pragma omp for nowait  #pragma omp for nowait
923                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
924                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
925                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
926                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
927                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 660  void Rectangle::setToIntegrals(vector<do Line 930  void Rectangle::setToIntegrals(vector<do
930    
931              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
932  #pragma omp for nowait  #pragma omp for nowait
933                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
934                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
935                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
936                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
937                      }  /* end of component loop i */                      }  /* end of component loop i */
# Line 672  void Rectangle::setToIntegrals(vector<do Line 942  void Rectangle::setToIntegrals(vector<do
942              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
943                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
944          } // end of parallel section          } // end of parallel section
     } else {  
         stringstream msg;  
         msg << "setToIntegrals() not implemented for "  
             << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
 }  
   
 void Rectangle::setToNormal(escript::Data& out) const  
 {  
     if (out.getFunctionSpace().getTypeCode() == FaceElements) {  
         out.requireWrite();  
 #pragma omp parallel  
         {  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);  
                     // set vector at two quadrature points  
                     *o++ = -1.;  
                     *o++ = 0.;  
                     *o++ = -1.;  
                     *o = 0.;  
                 }  
             }  
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);  
                     // set vector at two quadrature points  
                     *o++ = 1.;  
                     *o++ = 0.;  
                     *o++ = 1.;  
                     *o = 0.;  
                 }  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);  
                     // set vector at two quadrature points  
                     *o++ = 0.;  
                     *o++ = -1.;  
                     *o++ = 0.;  
                     *o = -1.;  
                 }  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);  
                     // set vector at two quadrature points  
                     *o++ = 0.;  
                     *o++ = 1.;  
                     *o++ = 0.;  
                     *o = 1.;  
                 }  
             }  
         } // end of parallel section  
     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
         out.requireWrite();  
 #pragma omp parallel  
         {  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);  
                     *o++ = -1.;  
                     *o = 0.;  
                 }  
             }  
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);  
                     *o++ = 1.;  
                     *o = 0.;  
                 }  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);  
                     *o++ = 0.;  
                     *o = -1.;  
                 }  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);  
                     *o++ = 0.;  
                     *o = 1.;  
                 }  
             }  
         } // end of parallel section  
   
     } else {  
         stringstream msg;  
         msg << "setToNormal() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
 }  
   
 void Rectangle::setToSize(escript::Data& out) const  
 {  
     if (out.getFunctionSpace().getTypeCode() == Elements  
             || out.getFunctionSpace().getTypeCode() == ReducedElements) {  
         out.requireWrite();  
         const dim_t numQuad=out.getNumDataPointsPerSample();  
         const double hSize=getFirstCoordAndSpacing(0).second;  
         const double vSize=getFirstCoordAndSpacing(1).second;  
         const double size=min(hSize,vSize);  
 #pragma omp parallel for  
         for (index_t k = 0; k < getNumElements(); ++k) {  
             double* o = out.getSampleDataRW(k);  
             fill(o, o+numQuad, size);  
         }  
     } else if (out.getFunctionSpace().getTypeCode() == FaceElements  
             || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
         out.requireWrite();  
         const dim_t numQuad=out.getNumDataPointsPerSample();  
         const double hSize=getFirstCoordAndSpacing(0).second;  
         const double vSize=getFirstCoordAndSpacing(1).second;  
 #pragma omp parallel  
         {  
             if (m_faceOffset[0] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     double* o = out.getSampleDataRW(m_faceOffset[0]+k1);  
                     fill(o, o+numQuad, vSize);  
                 }  
             }  
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     double* o = out.getSampleDataRW(m_faceOffset[1]+k1);  
                     fill(o, o+numQuad, vSize);  
                 }  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[2]+k0);  
                     fill(o, o+numQuad, hSize);  
                 }  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                     double* o = out.getSampleDataRW(m_faceOffset[3]+k0);  
                     fill(o, o+numQuad, hSize);  
                 }  
             }  
         } // end of parallel section  
   
     } else {  
         stringstream msg;  
         msg << "setToSize() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
 }  
   
 Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  
                                                 bool reducedColOrder) const  
 {  
     if (reducedRowOrder || reducedColOrder)  
         throw RipleyException("getPattern() not implemented for reduced order");  
   
     return m_pattern;  
 }  
   
 void Rectangle::Print_Mesh_Info(const bool full) const  
 {  
     RipleyDomain::Print_Mesh_Info(full);  
     if (full) {  
         cout << "     Id  Coordinates" << endl;  
         cout.precision(15);  
         cout.setf(ios::scientific, ios::floatfield);  
         pair<double,double> xdx = getFirstCoordAndSpacing(0);  
         pair<double,double> ydy = getFirstCoordAndSpacing(1);  
         for (index_t i=0; i < getNumNodes(); i++) {  
             cout << "  " << setw(5) << m_nodeId[i]  
                 << "  " << xdx.first+(i%m_N0)*xdx.second  
                 << "  " << ydy.first+(i/m_N0)*ydy.second << endl;  
         }  
     }  
 }  
   
 IndexVector Rectangle::getNumNodesPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_N0);  
     ret.push_back(m_N1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumElementsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NE0);  
     ret.push_back(m_NE1);  
     return ret;  
 }  
   
 IndexVector Rectangle::getNumFacesPerBoundary() const  
 {  
     IndexVector ret(4, 0);  
     //left  
     if (m_offset0==0)  
         ret[0]=m_NE1;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         ret[1]=m_NE1;  
     //bottom  
     if (m_offset1==0)  
         ret[2]=m_NE0;  
     //top  
     if (m_mpiInfo->rank/m_NX==m_NY-1)  
         ret[3]=m_NE0;  
     return ret;  
 }  
   
 pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const  
 {  
     if (dim==0) {  
         return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);  
     } else if (dim==1) {  
         return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);  
     }  
     throw RipleyException("getFirstCoordAndSpacing(): invalid argument");  
 }  
   
 //protected  
 dim_t Rectangle::getNumDOF() const  
 {  
     return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;  
 }  
   
 //protected  
 dim_t Rectangle::getNumFaceElements() const  
 {  
     const IndexVector faces = getNumFacesPerBoundary();  
     dim_t n=0;  
     for (size_t i=0; i<faces.size(); i++)  
         n+=faces[i];  
     return n;  
 }  
   
 //protected  
 void Rectangle::assembleCoordinates(escript::Data& arg) const  
 {  
     escriptDataC x = arg.getDataC();  
     int numDim = m_numDim;  
     if (!isDataPointShapeEqual(&x, 1, &numDim))  
         throw RipleyException("setToX: Invalid Data object shape");  
     if (!numSamplesEqual(&x, 1, getNumNodes()))  
         throw RipleyException("setToX: Illegal number of samples in Data object");  
   
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
     arg.requireWrite();  
 #pragma omp parallel for  
     for (dim_t i1 = 0; i1 < m_N1; i1++) {  
         for (dim_t i0 = 0; i0 < m_N0; i0++) {  
             double* point = arg.getSampleDataRW(i0+m_N0*i1);  
             point[0] = xdx.first+i0*xdx.second;  
             point[1] = ydy.first+i1*ydy.second;  
         }  
945      }      }
946  }  }
947    
# Line 1279  void Rectangle::interpolateNodesOnElemen Line 1269  void Rectangle::interpolateNodesOnElemen
1269      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1270      if (reduced) {      if (reduced) {
1271          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */  
1272          const double c0 = .25;          const double c0 = .25;
1273  #pragma omp parallel for  #pragma omp parallel for
1274          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 1294  void Rectangle::interpolateNodesOnElemen Line 1283  void Rectangle::interpolateNodesOnElemen
1283                  } /* end of component loop i */                  } /* end of component loop i */
1284              } /* end of k0 loop */              } /* end of k0 loop */
1285          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */  
1286      } else {      } else {
1287          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */  
1288          const double c0 = .16666666666666666667;          const double c0 = .16666666666666666667;
1289          const double c1 = .044658198738520451079;          const double c1 = .044658198738520451079;
1290          const double c2 = .62200846792814621559;          const double c2 = .62200846792814621559;
# Line 1317  void Rectangle::interpolateNodesOnElemen Line 1304  void Rectangle::interpolateNodesOnElemen
1304                  } /* end of component loop i */                  } /* end of component loop i */
1305              } /* end of k0 loop */              } /* end of k0 loop */
1306          } /* end of k1 loop */          } /* end of k1 loop */
         /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */  
1307      }      }
1308  }  }
1309    
# Line 1331  void Rectangle::interpolateNodesOnFaces( Line 1317  void Rectangle::interpolateNodesOnFaces(
1317          const double c0 = .5;          const double c0 = .5;
1318  #pragma omp parallel  #pragma omp parallel
1319          {          {
             /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */  
1320              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1321  #pragma omp for nowait  #pragma omp for nowait
1322                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 1376  void Rectangle::interpolateNodesOnFaces( Line 1361  void Rectangle::interpolateNodesOnFaces(
1361                      } /* end of component loop i */                      } /* end of component loop i */
1362                  } /* end of k0 loop */                  } /* end of k0 loop */
1363              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */  
1364          } // end of parallel section          } // end of parallel section
1365      } else {      } else {
1366          out.requireWrite();          out.requireWrite();
# Line 1384  void Rectangle::interpolateNodesOnFaces( Line 1368  void Rectangle::interpolateNodesOnFaces(
1368          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1369  #pragma omp parallel  #pragma omp parallel
1370          {          {
             /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */  
1371              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1372  #pragma omp for nowait  #pragma omp for nowait
1373                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 1433  void Rectangle::interpolateNodesOnFaces( Line 1416  void Rectangle::interpolateNodesOnFaces(
1416                      } /* end of component loop i */                      } /* end of component loop i */
1417                  } /* end of k0 loop */                  } /* end of k0 loop */
1418              } /* end of face 3 */              } /* end of face 3 */
             /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */  
1419          } // end of parallel section          } // end of parallel section
1420      }      }
1421  }  }
# Line 1447  void Rectangle::assemblePDESingle(Paso_S Line 1429  void Rectangle::assemblePDESingle(Paso_S
1429  {  {
1430      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1431      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
     /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */  
1432      const double w0 = -0.1555021169820365539*h1/h0;      const double w0 = -0.1555021169820365539*h1/h0;
1433      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
1434        const double w2 = -0.15550211698203655390;
1435        const double w3 = 0.041666666666666666667*h0/h1;
1436        const double w4 = 0.15550211698203655390;
1437        const double w5 = -0.041666666666666666667;
1438        const double w6 = -0.01116454968463011277*h1/h0;
1439        const double w7 = 0.011164549684630112770;
1440        const double w8 = -0.011164549684630112770;
1441        const double w9 = -0.041666666666666666667*h1/h0;
1442      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.041666666666666666667*h0/h1;
1443      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = 0.1555021169820365539*h1/h0;
1444      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = 0.1555021169820365539*h0/h1;
# Line 1460  void Rectangle::assemblePDESingle(Paso_S Line 1449  void Rectangle::assemblePDESingle(Paso_S
1449      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*h0/h1;
1450      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*h1/h0;
1451      const double w19 = 0.25000000000000000000;      const double w19 = 0.25000000000000000000;
     const double w2 = -0.15550211698203655390;  
1452      const double w20 = -0.25000000000000000000;      const double w20 = -0.25000000000000000000;
1453      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*h0/h1;
1454      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*h1/h0;
# Line 1471  void Rectangle::assemblePDESingle(Paso_S Line 1459  void Rectangle::assemblePDESingle(Paso_S
1459      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = -0.33333333333333333333*h0/h1;
1460      const double w28 = -0.032861463941450536761*h1;      const double w28 = -0.032861463941450536761*h1;
1461      const double w29 = -0.032861463941450536761*h0;      const double w29 = -0.032861463941450536761*h0;
     const double w3 = 0.041666666666666666667*h0/h1;  
1462      const double w30 = -0.12264065304058601714*h1;      const double w30 = -0.12264065304058601714*h1;
1463      const double w31 = -0.0023593469594139828636*h1;      const double w31 = -0.0023593469594139828636*h1;
1464      const double w32 = -0.008805202725216129906*h0;      const double w32 = -0.008805202725216129906*h0;
# Line 1482  void Rectangle::assemblePDESingle(Paso_S Line 1469  void Rectangle::assemblePDESingle(Paso_S
1469      const double w37 = 0.0023593469594139828636*h1;      const double w37 = 0.0023593469594139828636*h1;
1470      const double w38 = 0.12264065304058601714*h1;      const double w38 = 0.12264065304058601714*h1;
1471      const double w39 = 0.032861463941450536761*h0;      const double w39 = 0.032861463941450536761*h0;
     const double w4 = 0.15550211698203655390;  
1472      const double w40 = -0.12264065304058601714*h0;      const double w40 = -0.12264065304058601714*h0;
1473      const double w41 = -0.0023593469594139828636*h0;      const double w41 = -0.0023593469594139828636*h0;
1474      const double w42 = 0.0023593469594139828636*h0;      const double w42 = 0.0023593469594139828636*h0;
# Line 1493  void Rectangle::assemblePDESingle(Paso_S Line 1479  void Rectangle::assemblePDESingle(Paso_S
1479      const double w47 = 0.16666666666666666667*h1;      const double w47 = 0.16666666666666666667*h1;
1480      const double w48 = 0.083333333333333333333*h0;      const double w48 = 0.083333333333333333333*h0;
1481      const double w49 = -0.16666666666666666667*h0;      const double w49 = -0.16666666666666666667*h0;
     const double w5 = -0.041666666666666666667;  
1482      const double w50 = 0.16666666666666666667*h0;      const double w50 = 0.16666666666666666667*h0;
1483      const double w51 = -0.083333333333333333333*h1;      const double w51 = -0.083333333333333333333*h1;
1484      const double w52 = 0.025917019497006092316*h0*h1;      const double w52 = 0.025917019497006092316*h0*h1;
# Line 1504  void Rectangle::assemblePDESingle(Paso_S Line 1489  void Rectangle::assemblePDESingle(Paso_S
1489      const double w57 = 0.055555555555555555556*h0*h1;      const double w57 = 0.055555555555555555556*h0*h1;
1490      const double w58 = 0.027777777777777777778*h0*h1;      const double w58 = 0.027777777777777777778*h0*h1;
1491      const double w59 = 0.11111111111111111111*h0*h1;      const double w59 = 0.11111111111111111111*h0*h1;
     const double w6 = -0.01116454968463011277*h1/h0;  
1492      const double w60 = -0.19716878364870322056*h1;      const double w60 = -0.19716878364870322056*h1;
1493      const double w61 = -0.19716878364870322056*h0;      const double w61 = -0.19716878364870322056*h0;
1494      const double w62 = -0.052831216351296779436*h0;      const double w62 = -0.052831216351296779436*h0;
# Line 1515  void Rectangle::assemblePDESingle(Paso_S Line 1499  void Rectangle::assemblePDESingle(Paso_S
1499      const double w67 = 0.052831216351296779436*h0;      const double w67 = 0.052831216351296779436*h0;
1500      const double w68 = -0.5*h1;      const double w68 = -0.5*h1;
1501      const double w69 = -0.5*h0;      const double w69 = -0.5*h0;
     const double w7 = 0.011164549684630112770;  
1502      const double w70 = 0.5*h1;      const double w70 = 0.5*h1;
1503      const double w71 = 0.5*h0;      const double w71 = 0.5*h0;
1504      const double w72 = 0.1555021169820365539*h0*h1;      const double w72 = 0.1555021169820365539*h0*h1;
1505      const double w73 = 0.041666666666666666667*h0*h1;      const double w73 = 0.041666666666666666667*h0*h1;
1506      const double w74 = 0.01116454968463011277*h0*h1;      const double w74 = 0.01116454968463011277*h0*h1;
1507      const double w75 = 0.25*h0*h1;      const double w75 = 0.25*h0*h1;
     const double w8 = -0.011164549684630112770;  
     const double w9 = -0.041666666666666666667*h1/h0;  
     /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */  
1508    
1509      rhs.requireWrite();      rhs.requireWrite();
1510  #pragma omp parallel  #pragma omp parallel
# Line 1538  void Rectangle::assemblePDESingle(Paso_S Line 1518  void Rectangle::assemblePDESingle(Paso_S
1518                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1519                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1520                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
                     /*** GENERATOR SNIP_PDE_SINGLE TOP */  
1521                      ///////////////                      ///////////////
1522                      // process A //                      // process A //
1523                      ///////////////                      ///////////////
# Line 1547  void Rectangle::assemblePDESingle(Paso_S Line 1526  void Rectangle::assemblePDESingle(Paso_S
1526                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1527                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
1528                              const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];                              const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
                             const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];  
1529                              const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];                              const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1530                                const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1531                              const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];                              const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1532                              const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];                              const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
                             const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];  
1533                              const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];                              const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1534                                const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1535                              const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];                              const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1536                              const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];                              const register double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
                             const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];  
1537                              const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];                              const register double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1538                                const register double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1539                              const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];                              const register double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1540                              const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];                              const register double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
                             const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];  
1541                              const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];                              const register double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1542                                const register double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1543                              const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];                              const register double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1544                                const register double tmp0_0 = A_01_0 + A_01_3;
1545                                const register double tmp1_0 = A_00_0 + A_00_1;
1546                                const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
1547                                const register double tmp3_0 = A_00_2 + A_00_3;
1548                              const register double tmp4_0 = A_10_1 + A_10_2;                              const register double tmp4_0 = A_10_1 + A_10_2;
1549                                const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
1550                                const register double tmp6_0 = A_01_3 + A_10_0;
1551                                const register double tmp7_0 = A_01_0 + A_10_3;
1552                                const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
1553                                const register double tmp9_0 = A_01_0 + A_10_0;
1554                              const register double tmp12_0 = A_11_0 + A_11_2;                              const register double tmp12_0 = A_11_0 + A_11_2;
                             const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;  
1555                              const register double tmp10_0 = A_01_3 + A_10_3;                              const register double tmp10_0 = A_01_3 + A_10_3;
1556                              const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                              const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
                             const register double tmp0_0 = A_01_0 + A_01_3;  
1557                              const register double tmp13_0 = A_01_2 + A_10_1;                              const register double tmp13_0 = A_01_2 + A_10_1;
                             const register double tmp3_0 = A_00_2 + A_00_3;  
1558                              const register double tmp11_0 = A_11_1 + A_11_3;                              const register double tmp11_0 = A_11_1 + A_11_3;
1559                              const register double tmp18_0 = A_01_1 + A_10_1;                              const register double tmp18_0 = A_01_1 + A_10_1;
                             const register double tmp1_0 = A_00_0 + A_00_1;  
1560                              const register double tmp15_0 = A_01_1 + A_10_2;                              const register double tmp15_0 = A_01_1 + A_10_2;
                             const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;  
1561                              const register double tmp16_0 = A_10_0 + A_10_3;                              const register double tmp16_0 = A_10_0 + A_10_3;
                             const register double tmp6_0 = A_01_3 + A_10_0;  
1562                              const register double tmp17_0 = A_01_1 + A_01_2;                              const register double tmp17_0 = A_01_1 + A_01_2;
                             const register double tmp9_0 = A_01_0 + A_10_0;  
                             const register double tmp7_0 = A_01_0 + A_10_3;  
                             const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;  
1563                              const register double tmp19_0 = A_01_2 + A_10_2;                              const register double tmp19_0 = A_01_2 + A_10_2;
1564                                const register double tmp0_1 = A_10_3*w8;
1565                                const register double tmp1_1 = tmp0_0*w1;
1566                                const register double tmp2_1 = A_01_1*w4;
1567                                const register double tmp3_1 = tmp1_0*w0;
1568                                const register double tmp4_1 = A_01_2*w7;
1569                                const register double tmp5_1 = tmp2_0*w3;
1570                                const register double tmp6_1 = tmp3_0*w6;
1571                                const register double tmp7_1 = A_10_0*w2;
1572                                const register double tmp8_1 = tmp4_0*w5;
1573                                const register double tmp9_1 = tmp2_0*w10;
1574                              const register double tmp14_1 = A_10_0*w8;                              const register double tmp14_1 = A_10_0*w8;
1575                              const register double tmp23_1 = tmp3_0*w14;                              const register double tmp23_1 = tmp3_0*w14;
1576                              const register double tmp35_1 = A_01_0*w8;                              const register double tmp35_1 = A_01_0*w8;
1577                              const register double tmp54_1 = tmp13_0*w8;                              const register double tmp54_1 = tmp13_0*w8;
1578                              const register double tmp20_1 = tmp9_0*w4;                              const register double tmp20_1 = tmp9_0*w4;
1579                              const register double tmp25_1 = tmp12_0*w12;                              const register double tmp25_1 = tmp12_0*w12;
                             const register double tmp2_1 = A_01_1*w4;  
1580                              const register double tmp44_1 = tmp7_0*w7;                              const register double tmp44_1 = tmp7_0*w7;
1581                              const register double tmp26_1 = tmp10_0*w4;                              const register double tmp26_1 = tmp10_0*w4;
1582                              const register double tmp52_1 = tmp18_0*w8;                              const register double tmp52_1 = tmp18_0*w8;
1583                              const register double tmp48_1 = A_10_1*w7;                              const register double tmp48_1 = A_10_1*w7;
1584                              const register double tmp46_1 = A_01_3*w8;                              const register double tmp46_1 = A_01_3*w8;
1585                              const register double tmp50_1 = A_01_0*w2;                              const register double tmp50_1 = A_01_0*w2;
                             const register double tmp8_1 = tmp4_0*w5;  
1586                              const register double tmp56_1 = tmp19_0*w8;                              const register double tmp56_1 = tmp19_0*w8;
                             const register double tmp9_1 = tmp2_0*w10;  
1587                              const register double tmp19_1 = A_10_3*w2;                              const register double tmp19_1 = A_10_3*w2;
1588                              const register double tmp47_1 = A_10_2*w4;                              const register double tmp47_1 = A_10_2*w4;
1589                              const register double tmp16_1 = tmp3_0*w0;                              const register double tmp16_1 = tmp3_0*w0;
# Line 1610  void Rectangle::assemblePDESingle(Paso_S Line 1596  void Rectangle::assemblePDESingle(Paso_S
1596                              const register double tmp34_1 = tmp15_0*w8;                              const register double tmp34_1 = tmp15_0*w8;
1597                              const register double tmp33_1 = tmp14_0*w5;                              const register double tmp33_1 = tmp14_0*w5;
1598                              const register double tmp24_1 = tmp11_0*w13;                              const register double tmp24_1 = tmp11_0*w13;
                             const register double tmp3_1 = tmp1_0*w0;  
                             const register double tmp5_1 = tmp2_0*w3;  
1599                              const register double tmp43_1 = tmp17_0*w5;                              const register double tmp43_1 = tmp17_0*w5;
1600                              const register double tmp15_1 = A_01_2*w4;                              const register double tmp15_1 = A_01_2*w4;
1601                              const register double tmp53_1 = tmp19_0*w2;                              const register double tmp53_1 = tmp19_0*w2;
# Line 1630  void Rectangle::assemblePDESingle(Paso_S Line 1614  void Rectangle::assemblePDESingle(Paso_S
1614                              const register double tmp42_1 = tmp12_0*w16;                              const register double tmp42_1 = tmp12_0*w16;
1615                              const register double tmp49_1 = tmp12_0*w17;                              const register double tmp49_1 = tmp12_0*w17;
1616                              const register double tmp21_1 = tmp1_0*w11;                              const register double tmp21_1 = tmp1_0*w11;
                             const register double tmp1_1 = tmp0_0*w1;  
1617                              const register double tmp45_1 = tmp6_0*w4;                              const register double tmp45_1 = tmp6_0*w4;
                             const register double tmp7_1 = A_10_0*w2;  
                             const register double tmp6_1 = tmp3_0*w6;  
1618                              const register double tmp13_1 = tmp8_0*w1;                              const register double tmp13_1 = tmp8_0*w1;
1619                              const register double tmp36_1 = tmp16_0*w1;                              const register double tmp36_1 = tmp16_0*w1;
1620                              const register double tmp41_1 = A_01_3*w2;                              const register double tmp41_1 = A_01_3*w2;
1621                              const register double tmp30_1 = tmp12_0*w13;                              const register double tmp30_1 = tmp12_0*w13;
                             const register double tmp4_1 = A_01_2*w7;  
                             const register double tmp0_1 = A_10_3*w8;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
1622                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;                              EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
1623                              EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
1624                                EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
1625                              EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;                              EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
1626                              EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
1627                                EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
1628                              EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;                              EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
1629                                EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
1630                              EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1631                              EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
                             EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
1632                              EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;                              EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
1633                              EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;                              EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
1634                              EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;                              EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
1635                              EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;                              EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
1636                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
1637                                EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1638                            } else { // constant data
1639                              const register double A_00 = A_p[INDEX2(0,0,2)];                              const register double A_00 = A_p[INDEX2(0,0,2)];
                             const register double A_01 = A_p[INDEX2(0,1,2)];  
1640                              const register double A_10 = A_p[INDEX2(1,0,2)];                              const register double A_10 = A_p[INDEX2(1,0,2)];
1641                                const register double A_01 = A_p[INDEX2(0,1,2)];
1642                              const register double A_11 = A_p[INDEX2(1,1,2)];                              const register double A_11 = A_p[INDEX2(1,1,2)];
1643                              const register double tmp0_0 = A_01 + A_10;                              const register double tmp0_0 = A_01 + A_10;
1644                              const register double tmp0_1 = A_00*w18;                              const register double tmp0_1 = A_00*w18;
1645                              const register double tmp10_1 = A_01*w20;                              const register double tmp1_1 = A_01*w19;
1646                              const register double tmp12_1 = A_00*w26;                              const register double tmp2_1 = A_10*w20;
1647                                const register double tmp3_1 = A_11*w21;
1648                              const register double tmp4_1 = A_00*w22;                              const register double tmp4_1 = A_00*w22;
1649                                const register double tmp5_1 = tmp0_0*w19;
1650                                const register double tmp6_1 = A_11*w23;
1651                                const register double tmp7_1 = A_11*w25;
1652                              const register double tmp8_1 = A_00*w24;                              const register double tmp8_1 = A_00*w24;
                             const register double tmp13_1 = A_10*w19;  
1653                              const register double tmp9_1 = tmp0_0*w20;                              const register double tmp9_1 = tmp0_0*w20;
1654                              const register double tmp3_1 = A_11*w21;                              const register double tmp10_1 = A_01*w20;
1655                              const register double tmp11_1 = A_11*w27;                              const register double tmp11_1 = A_11*w27;
1656                              const register double tmp1_1 = A_01*w19;                              const register double tmp12_1 = A_00*w26;
1657                              const register double tmp6_1 = A_11*w23;                              const register double tmp13_1 = A_10*w19;
                             const register double tmp7_1 = A_11*w25;  
                             const register double tmp2_1 = A_10*w20;  
                             const register double tmp5_1 = tmp0_0*w19;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
1658                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1659                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1660                                EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1661                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1662                              EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1663                                EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1664                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
1665                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1666                              EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;                              EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
1667                              EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
                             EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
1668                              EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
1669                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                              EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
1670                              EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
1671                              EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
1672                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
1673                                EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
1674                          }                          }
1675                      }                      }
1676                      ///////////////                      ///////////////
# Line 1709  void Rectangle::assemblePDESingle(Paso_S Line 1688  void Rectangle::assemblePDESingle(Paso_S
1688                              const register double B_1_2 = B_p[INDEX2(1,2,2)];                              const register double B_1_2 = B_p[INDEX2(1,2,2)];
1689                              const register double B_0_3 = B_p[INDEX2(0,3,2)];                              const register double B_0_3 = B_p[INDEX2(0,3,2)];
1690                              const register double B_1_3 = B_p[INDEX2(1,3,2)];                              const register double B_1_3 = B_p[INDEX2(1,3,2)];
1691                              const register double tmp3_0 = B_0_0 + B_0_2;                              const register double tmp0_0 = B_1_0 + B_1_1;
1692                              const register double tmp1_0 = B_1_2 + B_1_3;                              const register double tmp1_0 = B_1_2 + B_1_3;
1693                              const register double tmp2_0 = B_0_1 + B_0_3;                              const register double tmp2_0 = B_0_1 + B_0_3;
1694                              const register double tmp0_0 = B_1_0 + B_1_1;                              const register double tmp3_0 = B_0_0 + B_0_2;
1695                              const register double tmp63_1 = B_1_1*w42;                              const register double tmp63_1 = B_1_1*w42;
1696                              const register double tmp79_1 = B_1_1*w40;                              const register double tmp79_1 = B_1_1*w40;
1697                              const register double tmp37_1 = tmp3_0*w35;                              const register double tmp37_1 = tmp3_0*w35;
# Line 1793  void Rectangle::assemblePDESingle(Paso_S Line 1772  void Rectangle::assemblePDESingle(Paso_S
1772                              const register double tmp78_1 = B_0_0*w34;                              const register double tmp78_1 = B_0_0*w34;
1773                              const register double tmp12_1 = tmp0_0*w36;                              const register double tmp12_1 = tmp0_0*w36;
1774                              const register double tmp75_1 = B_1_0*w32;                              const register double tmp75_1 = B_1_0*w32;
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
1775                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1776                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1777                                EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1778                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;
1779                              EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1780                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1781                              EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;                              EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1782                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1783                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1784                              EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
1785                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1786                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                              EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1787                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;
1788                              EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;                              EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1789                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1790                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1791                            } else { // constant data
1792                              const register double B_0 = B_p[0];                              const register double B_0 = B_p[0];
1793                              const register double B_1 = B_p[1];                              const register double B_1 = B_p[1];
                             const register double tmp6_1 = B_1*w50;  
                             const register double tmp1_1 = B_1*w45;  
                             const register double tmp5_1 = B_1*w49;  
                             const register double tmp4_1 = B_1*w48;  
1794                              const register double tmp0_1 = B_0*w44;                              const register double tmp0_1 = B_0*w44;
1795                                const register double tmp1_1 = B_1*w45;
1796                              const register double tmp2_1 = B_0*w46;                              const register double tmp2_1 = B_0*w46;
                             const register double tmp7_1 = B_0*w51;  
1797                              const register double tmp3_1 = B_0*w47;                              const register double tmp3_1 = B_0*w47;
1798                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const register double tmp4_1 = B_1*w48;
1799                              EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;                              const register double tmp5_1 = B_1*w49;
1800                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;                              const register double tmp6_1 = B_1*w50;
1801                                const register double tmp7_1 = B_0*w51;
1802                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;
1803                              EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
1804                                EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;
1805                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
1806                              EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1807                                EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
1808                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
1809                                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;
1810                              EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;                              EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;
1811                              EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;                              EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
                             EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;  
1812                              EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;                              EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;
1813                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
1814                              EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;                              EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;
1815                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;
1816                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;
1817                                EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;
1818                          }                          }
1819                      }                      }
1820                      ///////////////                      ///////////////
# Line 1853  void Rectangle::assemblePDESingle(Paso_S Line 1832  void Rectangle::assemblePDESingle(Paso_S
1832                              const register double C_1_2 = C_p[INDEX2(1,2,2)];                              const register double C_1_2 = C_p[INDEX2(1,2,2)];
1833                              const register double C_0_3 = C_p[INDEX2(0,3,2)];                              const register double C_0_3 = C_p[INDEX2(0,3,2)];
1834                              const register double C_1_3 = C_p[INDEX2(1,3,2)];                              const register double C_1_3 = C_p[INDEX2(1,3,2)];
1835                              const register double tmp2_0 = C_0_1 + C_0_3;                              const register double tmp0_0 = C_1_0 + C_1_1;
1836                              const register double tmp1_0 = C_1_2 + C_1_3;                              const register double tmp1_0 = C_1_2 + C_1_3;
1837                                const register double tmp2_0 = C_0_1 + C_0_3;
1838                              const register double tmp3_0 = C_0_0 + C_0_2;                              const register double tmp3_0 = C_0_0 + C_0_2;
                             const register double tmp0_0 = C_1_0 + C_1_1;  
1839                              const register double tmp64_1 = C_0_2*w30;                              const register double tmp64_1 = C_0_2*w30;
1840                              const register double tmp14_1 = C_0_2*w28;                              const register double tmp14_1 = C_0_2*w28;
1841                              const register double tmp19_1 = C_0_3*w31;                              const register double tmp19_1 = C_0_3*w31;
# Line 1937  void Rectangle::assemblePDESingle(Paso_S Line 1916  void Rectangle::assemblePDESingle(Paso_S
1916                              const register double tmp78_1 = C_0_0*w34;                              const register double tmp78_1 = C_0_0*w34;
1917                              const register double tmp11_1 = tmp0_0*w36;                              const register double tmp11_1 = tmp0_0*w36;
1918                              const register double tmp23_1 = C_0_2*w33;                              const register double tmp23_1 = C_0_2*w33;
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
1919                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                              EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
1920                              EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;
1921                                EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
1922                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;                              EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;
1923                              EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1924                                EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
1925                              EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;                              EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;
1926                                EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1927                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;                              EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;
1928                              EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                              EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
1929                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;                              EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;
1930                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;                              EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
1931                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;                              EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;
1932                              EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;                              EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;
1933                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
1934                                EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
1935                            } else { // constant data
1936                              const register double C_0 = C_p[0];                              const register double C_0 = C_p[0];
1937                              const register double C_1 = C_p[1];                              const register double C_1 = C_p[1];
1938                                const register double tmp0_1 = C_0*w47;
1939                              const register double tmp1_1 = C_1*w45;                              const register double tmp1_1 = C_1*w45;
1940                                const register double tmp2_1 = C_1*w48;
1941                              const register double tmp3_1 = C_0*w51;                              const register double tmp3_1 = C_0*w51;
1942                              const register double tmp4_1 = C_0*w44;                              const register double tmp4_1 = C_0*w44;
                             const register double tmp7_1 = C_0*w46;  
1943                              const register double tmp5_1 = C_1*w49;                              const register double tmp5_1 = C_1*w49;
                             const register double tmp2_1 = C_1*w48;  
                             const register double tmp0_1 = C_0*w47;  
1944                              const register double tmp6_1 = C_1*w50;                              const register double tmp6_1 = C_1*w50;
1945                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const register double tmp7_1 = C_0*w46;
                             EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;  
1946                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;
1947                              EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
1948                                EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;
1949                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
1950                              EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
1951                                EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;
1952                              EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;                              EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;
1953                                EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;
1954                              EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;                              EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;
1955                              EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
                             EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
1956                              EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;                              EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;
1957                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;                              EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;
1958                              EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;                              EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;
1959                              EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;
1960                                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
1961                                EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;
1962                          }                          }
1963                      }                      }
1964                      ///////////////                      ///////////////
# Line 1993  void Rectangle::assemblePDESingle(Paso_S Line 1972  void Rectangle::assemblePDESingle(Paso_S
1972                              const register double D_1 = D_p[1];                              const register double D_1 = D_p[1];
1973                              const register double D_2 = D_p[2];                              const register double D_2 = D_p[2];
1974                              const register double D_3 = D_p[3];                              const register double D_3 = D_p[3];
                             const register double tmp4_0 = D_1 + D_3;  
                             const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;  
                             const register double tmp5_0 = D_0 + D_2;  
1975                              const register double tmp0_0 = D_0 + D_1;                              const register double tmp0_0 = D_0 + D_1;
                             const register double tmp6_0 = D_0 + D_3;  
1976                              const register double tmp1_0 = D_2 + D_3;                              const register double tmp1_0 = D_2 + D_3;
1977                                const register double tmp2_0 = D_0 + D_1 + D_2 + D_3;
1978                              const register double tmp3_0 = D_1 + D_2;                              const register double tmp3_0 = D_1 + D_2;
1979                              const register double tmp16_1 = D_1*w56;                              const register double tmp4_0 = D_1 + D_3;
1980                              const register double tmp14_1 = tmp6_0*w54;                              const register double tmp5_0 = D_0 + D_2;
1981                              const register double tmp8_1 = D_3*w55;                              const register double tmp6_0 = D_0 + D_3;
1982                                const register double tmp0_1 = tmp0_0*w52;
1983                                const register double tmp1_1 = tmp1_0*w53;
1984                              const register double tmp2_1 = tmp2_0*w54;                              const register double tmp2_1 = tmp2_0*w54;
                             const register double tmp12_1 = tmp5_0*w52;  
                             const register double tmp4_1 = tmp0_0*w53;  
1985                              const register double tmp3_1 = tmp1_0*w52;                              const register double tmp3_1 = tmp1_0*w52;
1986                              const register double tmp13_1 = tmp4_0*w53;                              const register double tmp4_1 = tmp0_0*w53;
                             const register double tmp10_1 = tmp4_0*w52;  
                             const register double tmp1_1 = tmp1_0*w53;  
                             const register double tmp7_1 = D_3*w56;  
                             const register double tmp0_1 = tmp0_0*w52;  
                             const register double tmp11_1 = tmp5_0*w53;  
                             const register double tmp9_1 = D_0*w56;  
1987                              const register double tmp5_1 = tmp3_0*w54;                              const register double tmp5_1 = tmp3_0*w54;
                             const register double tmp18_1 = D_2*w56;  
                             const register double tmp17_1 = D_1*w55;  
1988                              const register double tmp6_1 = D_0*w55;                              const register double tmp6_1 = D_0*w55;
1989                                const register double tmp7_1 = D_3*w56;
1990                                const register double tmp8_1 = D_3*w55;
1991                                const register double tmp9_1 = D_0*w56;
1992                                const register double tmp10_1 = tmp4_0*w52;
1993                                const register double tmp11_1 = tmp5_0*w53;
1994                                const register double tmp12_1 = tmp5_0*w52;
1995                                const register double tmp13_1 = tmp4_0*w53;
1996                                const register double tmp14_1 = tmp6_0*w54;
1997                              const register double tmp15_1 = D_2*w55;                              const register double tmp15_1 = D_2*w55;
1998                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              const register double tmp16_1 = D_1*w56;
1999                              EM_S[INDEX2(1,2,4)]+=tmp2_1;                              const register double tmp17_1 = D_1*w55;
2000                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;                              const register double tmp18_1 = D_2*w56;
2001                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
2002                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;
2003                                EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;
2004                              EM_S[INDEX2(3,0,4)]+=tmp2_1;                              EM_S[INDEX2(3,0,4)]+=tmp2_1;
2005                              EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2006                                EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
2007                              EM_S[INDEX2(2,1,4)]+=tmp2_1;                              EM_S[INDEX2(2,1,4)]+=tmp2_1;
2008                                EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;
2009                              EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;                              EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;
2010                              EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;                              EM_S[INDEX2(1,2,4)]+=tmp2_1;
                             EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;  
2011                              EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;                              EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;
2012                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;
2013                              EM_S[INDEX2(0,3,4)]+=tmp2_1;                              EM_S[INDEX2(0,3,4)]+=tmp2_1;
2014                              EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;                              EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
2015                          } else { /* constant data */                              EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;
2016                              const register double D_0 = D_p[0];                              EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;
2017                              const register double tmp2_1 = D_0*w59;                          } else { // constant data
2018                              const register double tmp1_1 = D_0*w58;                              const register double tmp0_1 = D_p[0]*w57;
2019                              const register double tmp0_1 = D_0*w57;                              const register double tmp1_1 = D_p[0]*w58;
2020                              EM_S[INDEX2(0,1,4)]+=tmp0_1;                              const register double tmp2_1 = D_p[0]*w59;
                             EM_S[INDEX2(1,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1;  
2021                              EM_S[INDEX2(0,0,4)]+=tmp2_1;                              EM_S[INDEX2(0,0,4)]+=tmp2_1;
2022                              EM_S[INDEX2(3,3,4)]+=tmp2_1;                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
2023                                EM_S[INDEX2(2,0,4)]+=tmp0_1;
2024                              EM_S[INDEX2(3,0,4)]+=tmp1_1;                              EM_S[INDEX2(3,0,4)]+=tmp1_1;
2025                              EM_S[INDEX2(3,1,4)]+=tmp0_1;                              EM_S[INDEX2(0,1,4)]+=tmp0_1;
2026                                EM_S[INDEX2(1,1,4)]+=tmp2_1;
2027                              EM_S[INDEX2(2,1,4)]+=tmp1_1;                              EM_S[INDEX2(2,1,4)]+=tmp1_1;
2028                                EM_S[INDEX2(3,1,4)]+=tmp0_1;
2029                              EM_S[INDEX2(0,2,4)]+=tmp0_1;                              EM_S[INDEX2(0,2,4)]+=tmp0_1;
2030                              EM_S[INDEX2(2,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,2,4)]+=tmp1_1;
                             EM_S[INDEX2(1,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1;  
2031                              EM_S[INDEX2(2,2,4)]+=tmp2_1;                              EM_S[INDEX2(2,2,4)]+=tmp2_1;
2032                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(3,2,4)]+=tmp0_1;
2033                              EM_S[INDEX2(0,3,4)]+=tmp1_1;                              EM_S[INDEX2(0,3,4)]+=tmp1_1;
2034                              EM_S[INDEX2(1,1,4)]+=tmp2_1;                              EM_S[INDEX2(1,3,4)]+=tmp0_1;
2035                                EM_S[INDEX2(2,3,4)]+=tmp0_1;
2036                                EM_S[INDEX2(3,3,4)]+=tmp2_1;
2037                          }                          }
2038                      }                      }
2039                      ///////////////                      ///////////////
# Line 2073  void Rectangle::assemblePDESingle(Paso_S Line 2051  void Rectangle::assemblePDESingle(Paso_S
2051                              const register double X_1_2 = X_p[INDEX2(1,2,2)];                              const register double X_1_2 = X_p[INDEX2(1,2,2)];
2052                              const register double X_0_3 = X_p[INDEX2(0,3,2)];                              const register double X_0_3 = X_p[INDEX2(0,3,2)];
2053                              const register double X_1_3 = X_p[INDEX2(1,3,2)];                              const register double X_1_3 = X_p[INDEX2(1,3,2)];
2054                                const register double tmp0_0 = X_0_2 + X_0_3;
2055                              const register double tmp1_0 = X_1_1 + X_1_3;                              const register double tmp1_0 = X_1_1 + X_1_3;
                             const register double tmp3_0 = X_0_0 + X_0_1;  
2056                              const register double tmp2_0 = X_1_0 + X_1_2;                              const register double tmp2_0 = X_1_0 + X_1_2;
2057                              const register double tmp0_0 = X_0_2 + X_0_3;                              const register double tmp3_0 = X_0_0 + X_0_1;
2058                              const register double tmp8_1 = tmp2_0*w66;                              const register double tmp0_1 = tmp0_0*w63;
2059                              const register double tmp5_1 = tmp3_0*w64;                              const register double tmp1_1 = tmp1_0*w62;
2060                              const register double tmp14_1 = tmp0_0*w64;                              const register double tmp2_1 = tmp2_0*w61;
2061                              const register double tmp3_1 = tmp3_0*w60;                              const register double tmp3_1 = tmp3_0*w60;
2062                                const register double tmp4_1 = tmp0_0*w65;
2063                                const register double tmp5_1 = tmp3_0*w64;
2064                                const register double tmp6_1 = tmp2_0*w62;
2065                                const register double tmp7_1 = tmp1_0*w61;
2066                                const register double tmp8_1 = tmp2_0*w66;
2067                              const register double tmp9_1 = tmp3_0*w63;                              const register double tmp9_1 = tmp3_0*w63;
                             const register double tmp13_1 = tmp3_0*w65;  
                             const register double tmp12_1 = tmp1_0*w66;  
2068                              const register double tmp10_1 = tmp0_0*w60;                              const register double tmp10_1 = tmp0_0*w60;
                             const register double tmp2_1 = tmp2_0*w61;  
                             const register double tmp6_1 = tmp2_0*w62;  
                             const register double tmp4_1 = tmp0_0*w65;  
2069                              const register double tmp11_1 = tmp1_0*w67;                              const register double tmp11_1 = tmp1_0*w67;
2070                              const register double tmp1_1 = tmp1_0*w62;                              const register double tmp12_1 = tmp1_0*w66;
2071                              const register double tmp7_1 = tmp1_0*w61;                              const register double tmp13_1 = tmp3_0*w65;
2072                              const register double tmp0_1 = tmp0_0*w63;                              const register double tmp14_1 = tmp0_0*w64;
2073                              const register double tmp15_1 = tmp2_0*w67;                              const register double tmp15_1 = tmp2_0*w67;
2074                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2075                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;                              EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
2076                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                              EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
2077                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                              EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
2078                          } else { /* constant data */                          } else { // constant data
2079                              const register double X_0 = X_p[0];                              const register double tmp0_1 = X_p[1]*w69;
2080                              const register double X_1 = X_p[1];                              const register double tmp1_1 = X_p[0]*w68;
2081                              const register double tmp3_1 = X_1*w71;                              const register double tmp2_1 = X_p[0]*w70;
2082                              const register double tmp2_1 = X_0*w70;                              const register double tmp3_1 = X_p[1]*w71;
                             const register double tmp1_1 = X_0*w68;  
                             const register double tmp0_1 = X_1*w69;  
2083                              EM_F[0]+=tmp0_1 + tmp1_1;                              EM_F[0]+=tmp0_1 + tmp1_1;
2084                              EM_F[1]+=tmp0_1 + tmp2_1;                              EM_F[1]+=tmp0_1 + tmp2_1;
2085                              EM_F[2]+=tmp1_1 + tmp3_1;                              EM_F[2]+=tmp1_1 + tmp3_1;
# Line 2123  void Rectangle::assemblePDESingle(Paso_S Line 2099  void Rectangle::assemblePDESingle(Paso_S
2099                              const register double Y_3 = Y_p[3];                              const register double Y_3 = Y_p[3];
2100                              const register double tmp0_0 = Y_1 + Y_2;                              const register double tmp0_0 = Y_1 + Y_2;
2101                              const register double tmp1_0 = Y_0 + Y_3;                              const register double tmp1_0 = Y_0 + Y_3;
2102                              const register double tmp9_1 = Y_0*w74;                              const register double tmp0_1 = Y_0*w72;
2103                                const register double tmp1_1 = tmp0_0*w73;
2104                                const register double tmp2_1 = Y_3*w74;
2105                                const register double tmp3_1 = Y_1*w72;
2106                              const register double tmp4_1 = tmp1_0*w73;                              const register double tmp4_1 = tmp1_0*w73;
2107                              const register double tmp5_1 = Y_2*w74;                              const register double tmp5_1 = Y_2*w74;
                             const register double tmp7_1 = Y_1*w74;  
2108                              const register double tmp6_1 = Y_2*w72;                              const register double tmp6_1 = Y_2*w72;
2109                              const register double tmp2_1 = Y_3*w74;                              const register double tmp7_1 = Y_1*w74;
2110                              const register double tmp8_1 = Y_3*w72;                              const register double tmp8_1 = Y_3*w72;
2111                              const register double tmp3_1 = Y_1*w72;                              const register double tmp9_1 = Y_0*w74;
                             const register double tmp0_1 = Y_0*w72;  
                             const register double tmp1_1 = tmp0_0*w73;  
2112                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
2113                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;                              EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
2114                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;                              EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;
2115                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;                              EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;
2116                          } else { /* constant data */                          } else { // constant data
2117                              const register double Y_0 = Y_p[0];                              const register double tmp0_1 = Y_p[0]*w75;
                             const register double tmp0_1 = Y_0*w75;  
2118                              EM_F[0]+=tmp0_1;                              EM_F[0]+=tmp0_1;
2119                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
2120                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
2121                              EM_F[3]+=tmp0_1;                              EM_F[3]+=tmp0_1;
2122                          }                          }
2123                      }                      }
                     /* GENERATOR SNIP_PDE_SINGLE BOTTOM */  
2124    
2125                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2126                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
# Line 2156  void Rectangle::assemblePDESingle(Paso_S Line 2130  void Rectangle::assemblePDESingle(Paso_S
2130                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2131                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2132                      if (add_EM_F) {                      if (add_EM_F) {
                         //cout << "-- AddtoRHS -- " << endl;  
2133                          double *F_p=rhs.getSampleDataRW(0);                          double *F_p=rhs.getSampleDataRW(0);
2134                          for (index_t i=0; i<rowIndex.size(); i++) {                          for (index_t i=0; i<rowIndex.size(); i++) {
2135                              if (rowIndex[i]<getNumDOF()) {                              if (rowIndex[i]<getNumDOF()) {
2136                                  F_p[rowIndex[i]]+=EM_F[i];                                  F_p[rowIndex[i]]+=EM_F[i];
                                 //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;  
2137                              }                              }
2138                          }                          }
                         //cout << "---"<<endl;  
2139                      }                      }
2140                      if (add_EM_S) {                      if (add_EM_S) {
                         //cout << "-- AddtoSystem -- " << endl;  
2141                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2142                      }                      }
2143    
# Line 2186  void Rectangle::assemblePDESingleReduced Line 2156  void Rectangle::assemblePDESingleReduced
2156  {  {
2157      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
2158      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
2159      /*** GENERATOR SNIP_PDE_SINGLE_REDUCED_PRE TOP */      const double w0 = -.25*h1/h0;
2160      const double w0 = -0.25*h1/h0;      const double w1 = .25;
2161      const double w1 = 0.25;      const double w2 = -.25;
2162      const double w2 = -0.25;      const double w3 = .25*h0/h1;
2163      const double w3 = 0.25*h0/h1;      const double w4 = -.25*h0/h1;
2164      const double w4 = -0.25*h0/h1;      const double w5 = .25*h1/h0;
2165      const double w5 = 0.25*h1/h0;      const double w6 = -.125*h1;
2166      const double w6 = -0.125*h1;      const double w7 = -.125*h0;
2167      const double w7 = -0.125*h0;      const double w8 = .125*h1;
2168      const double w8 = 0.125*h1;      const double w9 = .125*h0;
2169      const double w9 = 0.125*h0;      const double w10 = .0625*h0*h1;
2170      const double w10 = 0.0625*h0*h1;      const double w11 = -.5*h1;
2171      const double w11 = -0.5*h1;      const double w12 = -.5*h0;
2172      const double w12 = -0.5*h0;      const double w13 = .5*h1;
2173      const double w13 = 0.5*h1;      const double w14 = .5*h0;
2174      const double w14 = 0.5*h0;      const double w15 = .25*h0*h1;
     const double w15 = 0.25*h0*h1;  
     /* GENERATOR SNIP_PDE_SINGLE_REDUCED_PRE BOTTOM */  
2175    
2176      rhs.requireWrite();      rhs.requireWrite();
2177  #pragma omp parallel  #pragma omp parallel
# Line 2217  void Rectangle::assemblePDESingleReduced Line 2185  void Rectangle::assemblePDESingleReduced
2185                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
2186                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
2187                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
                     /*** GENERATOR SNIP_PDE_SINGLE_REDUCED TOP */  
2188                      ///////////////                      ///////////////
2189                      // process A //                      // process A //
2190                      ///////////////                      ///////////////
# Line 2225  void Rectangle::assemblePDESingleReduced Line 2192  void Rectangle::assemblePDESingleReduced
2192                          add_EM_S=true;                          add_EM_S=true;
2193                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2194                          const register double A_00 = A_p[INDEX2(0,0,2)];                          const register double A_00 = A_p[INDEX2(0,0,2)];
                         const register double A_01 = A_p[INDEX2(0,1,2)];  
2195                          const register double A_10 = A_p[INDEX2(1,0,2)];                          const register double A_10 = A_p[INDEX2(1,0,2)];
2196                            const register double A_01 = A_p[INDEX2(0,1,2)];
2197                          const register double A_11 = A_p[INDEX2(1,1,2)];                          const register double A_11 = A_p[INDEX2(1,1,2)];
2198                          const register double tmp0_0 = A_01 + A_10;                          const register double tmp0_0 = A_01 + A_10;
2199                            const register double tmp0_1 = A_11*w3;
2200                            const register double tmp1_1 = A_00*w0;
2201                          const register double tmp2_1 = A_01*w1;                          const register double tmp2_1 = A_01*w1;
                         const register double tmp7_1 = tmp0_0*w2;  
2202                          const register double tmp3_1 = A_10*w2;                          const register double tmp3_1 = A_10*w2;
                         const register double tmp6_1 = A_00*w5;  
                         const register double tmp1_1 = A_00*w0;  
2203                          const register double tmp4_1 = tmp0_0*w1;                          const register double tmp4_1 = tmp0_0*w1;
                         const register double tmp9_1 = A_01*w2;  
2204                          const register double tmp5_1 = A_11*w4;                          const register double tmp5_1 = A_11*w4;
2205                            const register double tmp6_1 = A_00*w5;
2206                            const register double tmp7_1 = tmp0_0*w2;
2207                          const register double tmp8_1 = A_10*w1;                          const register double tmp8_1 = A_10*w1;
2208                          const register double tmp0_1 = A_11*w3;                          const register double tmp9_1 = A_01*w2;
                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
2209                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2210                          EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2211                            EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2212                          EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;                          EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2213                          EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2214                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2215                          EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;                          EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2216                            EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2217                          EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2218                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                          EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
                         EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
2219                          EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;                          EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2220                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2221                          EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;                          EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2222                          EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;                          EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2223                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2224                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2225                      }                      }
2226                      ///////////////                      ///////////////
2227                      // process B //                      // process B //
# Line 2262  void Rectangle::assemblePDESingleReduced Line 2229  void Rectangle::assemblePDESingleReduced
2229                      if (!B.isEmpty()) {                      if (!B.isEmpty()) {
2230                          add_EM_S=true;                          add_EM_S=true;
2231                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2232                          const register double B_0 = B_p[0];                          const register double tmp2_1 = B_p[0]*w8;
2233                          const register double B_1 = B_p[1];                          const register double tmp0_1 = B_p[0]*w6;
2234                          const register double tmp2_1 = B_0*w8;                          const register double tmp3_1 = B_p[1]*w9;
2235                          const register double tmp0_1 = B_0*w6;                          const register double tmp1_1 = B_p[1]*w7;
                         const register double tmp3_1 = B_1*w9;  
                         const register double tmp1_1 = B_1*w7;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;  
2236                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2237                          EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2238                            EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
2239                          EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;
2240                          EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2241                            EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;
2242                          EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;
2243                            EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;
2244                          EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;
2245                          EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
                         EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;  
2246                          EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;
2247                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2248                          EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;
2249                          EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;
2250                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;
2251                            EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
2252                      }                      }
2253                      ///////////////                      ///////////////
2254                      // process C //                      // process C //
# Line 2291  void Rectangle::assemblePDESingleReduced Line 2256  void Rectangle::assemblePDESingleReduced
2256                      if (!C.isEmpty()) {                      if (!C.isEmpty()) {
2257                          add_EM_S=true;                          add_EM_S=true;
2258                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2259                          const register double C_0 = C_p[0];                          const register double tmp1_1 = C_p[1]*w7;
2260                          const register double C_1 = C_p[1];                          const register double tmp0_1 = C_p[0]*w8;
2261                          const register double tmp1_1 = C_1*w7;                          const register double tmp3_1 = C_p[0]*w6;
2262                          const register double tmp0_1 = C_0*w8;                          const register double tmp2_1 = C_p[1]*w9;
                         const register double tmp3_1 = C_0*w6;  
                         const register double tmp2_1 = C_1*w9;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;  
2263                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2264                          EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2265                            EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
2266                          EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
2267                          EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2268                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;
2269                          EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;
2270                            EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;
2271                          EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;
2272                          EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
                         EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
2273                          EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
2274                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2275                          EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;                          EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;
2276                          EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;
2277                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2278                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;
2279                      }                      }
2280                      ///////////////                      ///////////////
2281                      // process D //                      // process D //
# Line 2320  void Rectangle::assemblePDESingleReduced Line 2283  void Rectangle::assemblePDESingleReduced
2283                      if (!D.isEmpty()) {                      if (!D.isEmpty()) {
2284                          add_EM_S=true;                          add_EM_S=true;
2285                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2286                          const register double D_0 = D_p[0];                          const register double tmp0_1 = D_p[0]*w10;
                         const register double tmp0_1 = D_0*w10;  
                         EM_S[INDEX2(0,1,4)]+=tmp0_1;  
                         EM_S[INDEX2(1,2,4)]+=tmp0_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp0_1;  
2287                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,0,4)]+=tmp0_1;
2288                          EM_S[INDEX2(3,3,4)]+=tmp0_1;                          EM_S[INDEX2(1,0,4)]+=tmp0_1;
2289                            EM_S[INDEX2(2,0,4)]+=tmp0_1;
2290                          EM_S[INDEX2(3,0,4)]+=tmp0_1;                          EM_S[INDEX2(3,0,4)]+=tmp0_1;
2291                          EM_S[INDEX2(3,1,4)]+=tmp0_1;                          EM_S[INDEX2(0,1,4)]+=tmp0_1;
2292                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
2293                          EM_S[INDEX2(2,1,4)]+=tmp0_1;                          EM_S[INDEX2(2,1,4)]+=tmp0_1;
2294                            EM_S[INDEX2(3,1,4)]+=tmp0_1;
2295                          EM_S[INDEX2(0,2,4)]+=tmp0_1;                          EM_S[INDEX2(0,2,4)]+=tmp0_1;
2296                          EM_S[INDEX2(2,0,4)]+=tmp0_1;                          EM_S[INDEX2(1,2,4)]+=tmp0_1;
                         EM_S[INDEX2(1,3,4)]+=tmp0_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1;  
2297                          EM_S[INDEX2(2,2,4)]+=tmp0_1;                          EM_S[INDEX2(2,2,4)]+=tmp0_1;
2298                          EM_S[INDEX2(1,0,4)]+=tmp0_1;                          EM_S[INDEX2(3,2,4)]+=tmp0_1;
2299                          EM_S[INDEX2(0,3,4)]+=tmp0_1;                          EM_S[INDEX2(0,3,4)]+=tmp0_1;
2300                          EM_S[INDEX2(1,1,4)]+=tmp0_1;                          EM_S[INDEX2(1,3,4)]+=tmp0_1;
2301                            EM_S[INDEX2(2,3,4)]+=tmp0_1;
2302                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
2303                      }                      }
2304                      ///////////////                      ///////////////
2305                      // process X //                      // process X //
# Line 2345  void Rectangle::assemblePDESingleReduced Line 2307  void Rectangle::assemblePDESingleReduced
2307                      if (!X.isEmpty()) {                      if (!X.isEmpty()) {
2308                          add_EM_F=true;                          add_EM_F=true;
2309                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2310                          const register double X_0 = X_p[0];                          const register double tmp0_1 = X_p[0]*w11;
2311                          const register double X_1 = X_p[1];                          const register double tmp2_1 = X_p[0]*w13;
2312                          const register double tmp0_1 = X_0*w11;                          const register double tmp1_1 = X_p[1]*w12;
2313                          const register double tmp2_1 = X_0*w13;                          const register double tmp3_1 = X_p[1]*w14;
                         const register double tmp1_1 = X_1*w12;  
                         const register double tmp3_1 = X_1*w14;  
2314                          EM_F[0]+=tmp0_1 + tmp1_1;                          EM_F[0]+=tmp0_1 + tmp1_1;
2315                          EM_F[1]+=tmp1_1 + tmp2_1;                          EM_F[1]+=tmp1_1 + tmp2_1;
2316                          EM_F[2]+=tmp0_1 + tmp3_1;                          EM_F[2]+=tmp0_1 + tmp3_1;
# Line 2362  void Rectangle::assemblePDESingleReduced Line 2322  void Rectangle::assemblePDESingleReduced
2322                      if (!Y.isEmpty()) {                      if (!Y.isEmpty()) {
2323                          add_EM_F=true;                          add_EM_F=true;
2324                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2325                          const register double Y_0 = Y_p[0];                          const register double tmp0_1 = Y_p[0]*w15;
                         const register double tmp0_1 = Y_0*w15;  
2326                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0_1;
2327                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0_1;
2328                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0_1;
2329                          EM_F[3]+=tmp0_1;                          EM_F[3]+=tmp0_1;
2330                      }                      }
                     /* GENERATOR SNIP_PDE_SINGLE_REDUCED BOTTOM */  
2331    
2332                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2333                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
# Line 2379  void Rectangle::assemblePDESingleReduced Line 2337  void Rectangle::assemblePDESingleReduced
2337                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2338                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2339                      if (add_EM_F) {                      if (add_EM_F) {
                         //cout << "-- AddtoRHS -- " << endl;  
2340                          double *F_p=rhs.getSampleDataRW(0);                          double *F_p=rhs.getSampleDataRW(0);
2341                          for (index_t i=0; i<rowIndex.size(); i++) {                          for (index_t i=0; i<rowIndex.size(); i++) {
2342                              if (rowIndex[i]<getNumDOF()) {                              if (rowIndex[i]<getNumDOF()) {
2343                                  F_p[rowIndex[i]]+=EM_F[i];                                  F_p[rowIndex[i]]+=EM_F[i];
                                 //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;  
2344                              }                              }
2345                          }                          }
                         //cout << "---"<<endl;  
2346                      }                      }
2347                      if (add_EM_S) {                      if (add_EM_S) {
                         //cout << "-- AddtoSystem -- " << endl;  
2348                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2349                      }                      }
2350    
# Line 2417  void Rectangle::assemblePDESystem(Paso_S Line 2371  void Rectangle::assemblePDESystem(Paso_S
2371          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
2372      }      }
2373    
     /*** GENERATOR SNIP_PDE_SYSTEM_PRE TOP */  
2374      const double w0 = -0.1555021169820365539*h1/h0;      const double w0 = -0.1555021169820365539*h1/h0;
2375      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
2376        const double w2 = -0.15550211698203655390;
2377        const double w3 = 0.041666666666666666667*h0/h1;
2378        const double w4 = 0.15550211698203655390;
2379        const double w5 = -0.041666666666666666667;
2380        const double w6 = -0.01116454968463011277*h1/h0;
2381        const double w7 = 0.011164549684630112770;
2382        const double w8 = -0.011164549684630112770;
2383        const double w9 = -0.041666666666666666667*h1/h0;
2384      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.041666666666666666667*h0/h1;
2385      const double w11 = 0.1555021169820365539*h1/h0;      const double w11 = 0.1555021169820365539*h1/h0;
2386      const double w12 = 0.1555021169820365539*h0/h1;      const double w12 = 0.1555021169820365539*h0/h1;
# Line 2430  void Rectangle::assemblePDESystem(Paso_S Line 2391  void Rectangle::assemblePDESystem(Paso_S
2391      const double w17 = -0.1555021169820365539*h0/h1;      const double w17 = -0.1555021169820365539*h0/h1;
2392      const double w18 = -0.33333333333333333333*h1/h0;      const double w18 = -0.33333333333333333333*h1/h0;
2393      const double w19 = 0.25000000000000000000;      const double w19 = 0.25000000000000000000;
     const double w2 = -0.15550211698203655390;  
2394      const double w20 = -0.25000000000000000000;      const double w20 = -0.25000000000000000000;
2395      const double w21 = 0.16666666666666666667*h0/h1;      const double w21 = 0.16666666666666666667*h0/h1;
2396      const double w22 = -0.16666666666666666667*h1/h0;      const double w22 = -0.16666666666666666667*h1/h0;
# Line 2441  void Rectangle::assemblePDESystem(Paso_S Line 2401  void Rectangle::assemblePDESystem(Paso_S
2401      const double w27 = -0.33333333333333333333*h0/h1;      const double w27 = -0.33333333333333333333*h0/h1;
2402      const double w28 = -0.032861463941450536761*h1;      const double w28 = -0.032861463941450536761*h1;
2403      const double w29 = -0.032861463941450536761*h0;      const double w29 = -0.032861463941450536761*h0;
     const double w3 = 0.041666666666666666667*h0/h1;  
2404      const double w30 = -0.12264065304058601714*h1;      const double w30 = -0.12264065304058601714*h1;
2405      const double w31 = -0.0023593469594139828636*h1;      const double w31 = -0.0023593469594139828636*h1;
2406      const double w32 = -0.008805202725216129906*h0;      const double w32 = -0.008805202725216129906*h0;
# Line 2452  void Rectangle::assemblePDESystem(Paso_S Line 2411  void Rectangle::assemblePDESystem(Paso_S
2411      const double w37 = 0.0023593469594139828636*h1;      const double w37 = 0.0023593469594139828636*h1;
2412      const double w38 = 0.12264065304058601714*h1;      const double w38 = 0.12264065304058601714*h1;
2413      const double w39 = 0.032861463941450536761*h0;      const double w39 = 0.032861463941450536761*h0;
     const double w4 = 0.15550211698203655390;  
2414      const double w40 = -0.12264065304058601714*h0;      const double w40 = -0.12264065304058601714*h0;
2415      const double w41 = -0.0023593469594139828636*h0;      const double w41 = -0.0023593469594139828636*h0;
2416      const double w42 = 0.0023593469594139828636*h0;      const double w42 = 0.0023593469594139828636*h0;
# Line 2463  void Rectangle::assemblePDESystem(Paso_S Line 2421  void Rectangle::assemblePDESystem(Paso_S
2421      const double w47 = 0.16666666666666666667*h1;      const double w47 = 0.16666666666666666667*h1;
2422      const double w48 = 0.083333333333333333333*h0;      const double w48 = 0.083333333333333333333*h0;
2423      const double w49 = -0.16666666666666666667*h0;      const double w49 = -0.16666666666666666667*h0;
     const double w5 = -0.041666666666666666667;  
2424      const double w50 = 0.16666666666666666667*h0;      const double w50 = 0.16666666666666666667*h0;
2425      const double w51 = -0.083333333333333333333*h1;      const double w51 = -0.083333333333333333333*h1;
2426      const double w52 = 0.025917019497006092316*h0*h1;      const double w52 = 0.025917019497006092316*h0*h1;
# Line 2474  void Rectangle::assemblePDESystem(Paso_S Line 2431  void Rectangle::assemblePDESystem(Paso_S
2431      const double w57 = 0.055555555555555555556*h0*h1;      const double w57 = 0.055555555555555555556*h0*h1;
2432      const double w58 = 0.027777777777777777778*h0*h1;      const double w58 = 0.027777777777777777778*h0*h1;
2433      const double w59 = 0.11111111111111111111*h0*h1;      const double w59 = 0.11111111111111111111*h0*h1;
     const double w6 = -0.01116454968463011277*h1/h0;  
2434      const double w60 = -0.19716878364870322056*h1;      const double w60 = -0.19716878364870322056*h1;
2435      const double w61 = -0.19716878364870322056*h0;      const double w61 = -0.19716878364870322056*h0;
2436      const double w62 = -0.052831216351296779436*h0;      const double w62 = -0.052831216351296779436*h0;
# Line 2485  void Rectangle::assemblePDESystem(Paso_S Line 2441  void Rectangle::assemblePDESystem(Paso_S
2441      const double w67 = 0.052831216351296779436*h0;      const double w67 = 0.052831216351296779436*h0;
2442      const double w68 = -0.5*h1;      const double w68 = -0.5*h1;
2443      const double w69 = -0.5*h0;      const double w69 = -0.5*h0;
     const double w7 = 0.011164549684630112770;  
2444      const double w70 = 0.5*h1;      const double w70 = 0.5*h1;
2445      const double w71 = 0.5*h0;      const double w71 = 0.5*h0;
2446      const double w72 = 0.1555021169820365539*h0*h1;      const double w72 = 0.1555021169820365539*h0*h1;
2447      const double w73 = 0.041666666666666666667*h0*h1;      const double w73 = 0.041666666666666666667*h0*h1;
2448      const double w74 = 0.01116454968463011277*h0*h1;      const double w74 = 0.01116454968463011277*h0*h1;
2449      const double w75 = 0.25*h0*h1;      const double w75 = 0.25*h0*h1;
     const double w8 = -0.011164549684630112770;  
     const double w9 = -0.041666666666666666667*h1/h0;  
     /* GENERATOR SNIP_PDE_SYSTEM_PRE BOTTOM */  
2450    
2451      rhs.requireWrite();      rhs.requireWrite();
2452  #pragma omp parallel  #pragma omp parallel
# Line 2508  void Rectangle::assemblePDESystem(Paso_S Line 2460  void Rectangle::assemblePDESystem(Paso_S
2460                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
2461                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
2462                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
                     /*** GENERATOR SNIP_PDE_SYSTEM TOP */  
2463                      ///////////////                      ///////////////
2464                      // process A //                      // process A //
2465                      ///////////////                      ///////////////
# Line 2534  void Rectangle::assemblePDESystem(Paso_S Line 2485  void Rectangle::assemblePDESystem(Paso_S
2485                                      const register double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];                                      const register double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];
2486                                      const register double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];                                      const register double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];
2487                                      const register double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];                                      const register double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];
                                     const register double tmp4_0 = A_10_1 + A_10_2;  
                                     const register double tmp12_0 = A_11_0 + A_11_2;  
                                     const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;  
                                     const register double tmp10_0 = A_01_3 + A_10_3;  
                                     const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;  
2488                                      const register double tmp0_0 = A_01_0 + A_01_3;                                      const register double tmp0_0 = A_01_0 + A_01_3;
                                     const register double tmp13_0 = A_01_2 + A_10_1;  
                                     const register double tmp3_0 = A_00_2 + A_00_3;  
                                     const register double tmp11_0 = A_11_1 + A_11_3;  
                                     const register double tmp18_0 = A_01_1 + A_10_1;  
2489                                      const register double tmp1_0 = A_00_0 + A_00_1;                                      const register double tmp1_0 = A_00_0 + A_00_1;
2490                                      const register double tmp15_0 = A_01_1 + A_10_2;                                      const register double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;
2491                                        const register double tmp3_0 = A_00_2 + A_00_3;
2492                                        const register double tmp4_0 = A_10_1 + A_10_2;
2493                                      const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                                      const register double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;
                                     const register double tmp16_0 = A_10_0 + A_10_3;  
2494                                      const register double tmp6_0 = A_01_3 + A_10_0;                                      const register double tmp6_0 = A_01_3 + A_10_0;
                                     const register double tmp17_0 = A_01_1 + A_01_2;  
                                     const register double tmp9_0 = A_01_0 + A_10_0;  
2495                                      const register double tmp7_0 = A_01_0 + A_10_3;                                      const register double tmp7_0 = A_01_0 + A_10_3;
2496                                      const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                                      const register double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;
2497                                        const register double tmp9_0 = A_01_0 + A_10_0;
2498                                        const register double tmp10_0 = A_01_3 + A_10_3;
2499                                        const register double tmp11_0 = A_11_1 + A_11_3;
2500                                        const register double tmp12_0 = A_11_0 + A_11_2;
2501                                        const register double tmp13_0 = A_01_2 + A_10_1;
2502                                        const register double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;
2503                                        const register double tmp15_0 = A_01_1 + A_10_2;
2504                                        const register double tmp16_0 = A_10_0 + A_10_3;
2505                                        const register double tmp17_0 = A_01_1 + A_01_2;
2506                                        const register double tmp18_0 = A_01_1 + A_10_1;
2507                                      const register double tmp19_0 = A_01_2 + A_10_2;                                      const register double tmp19_0 = A_01_2 + A_10_2;
2508                                      const register double tmp14_1 = A_10_0*w8;                                      const register double tmp0_1 = A_10_3*w8;
2509                                      const register double tmp23_1 = tmp3_0*w14;                                      const register double tmp1_1 = tmp0_0*w1;
                                     const register double tmp35_1 = A_01_0*w8;  
                                     const register double tmp54_1 = tmp13_0*w8;  
                                     const register double tmp20_1 = tmp9_0*w4;  
                                     const register double tmp25_1 = tmp12_0*w12;  
2510                                      const register double tmp2_1 = A_01_1*w4;                                      const register double tmp2_1 = A_01_1*w4;
2511                                      const register double tmp44_1 = tmp7_0*w7;                                      const register double tmp3_1 = tmp1_0*w0;
2512                                      const register double tmp26_1 = tmp10_0*w4;                                      const register double tmp4_1 = A_01_2*w7;
2513                                      const register double tmp52_1 = tmp18_0*w8;                                      const register double tmp5_1 = tmp2_0*w3;
2514                                      const register double tmp48_1 = A_10_1*w7;                                      const register double tmp6_1 = tmp3_0*w6;
2515                                      const register double tmp46_1 = A_01_3*w8;                                      const register double tmp7_1 = A_10_0*w2;
                                     const register double tmp50_1 = A_01_0*w2;  
2516                                      const register double tmp8_1 = tmp4_0*w5;                                      const register double tmp8_1 = tmp4_0*w5;
                                     const register double tmp56_1 = tmp19_0*w8;  
2517                                      const register double tmp9_1 = tmp2_0*w10;                                      const register double tmp9_1 = tmp2_0*w10;
2518                                      const register double tmp19_1 = A_10_3*w2;                                      const register double tmp10_1 = tmp5_0*w9;
2519                                      const register double tmp47_1 = A_10_2*w4;                                      const register double tmp11_1 = tmp6_0*w7;
2520                                        const register double tmp12_1 = tmp7_0*w4;
2521                                        const register double tmp13_1 = tmp8_0*w1;
2522                                        const register double tmp14_1 = A_10_0*w8;
2523                                        const register double tmp15_1 = A_01_2*w4;
2524                                      const register double tmp16_1 = tmp3_0*w0;                                      const register double tmp16_1 = tmp3_0*w0;
2525                                        const register double tmp17_1 = A_01_1*w7;
2526                                      const register double tmp18_1 = tmp1_0*w6;                                      const register double tmp18_1 = tmp1_0*w6;
2527                                      const register double tmp31_1 = tmp11_0*w12;                                      const register double tmp19_1 = A_10_3*w2;
2528                                      const register double tmp55_1 = tmp15_0*w2;                                      const register double tmp20_1 = tmp9_0*w4;
2529                                      const register double tmp39_1 = A_10_2*w7;                                      const register double tmp21_1 = tmp1_0*w11;
2530                                      const register double tmp11_1 = tmp6_0*w7;                                      const register double tmp22_1 = tmp10_0*w7;
2531                                      const register double tmp40_1 = tmp11_0*w17;                                      const register double tmp23_1 = tmp3_0*w14;
                                     const register double tmp34_1 = tmp15_0*w8;  
                                     const register double tmp33_1 = tmp14_0*w5;  
2532                                      const register double tmp24_1 = tmp11_0*w13;                                      const register double tmp24_1 = tmp11_0*w13;
2533                                      const register double tmp3_1 = tmp1_0*w0;                                      const register double tmp25_1 = tmp12_0*w12;
2534                                      const register double tmp5_1 = tmp2_0*w3;                                      const register double tmp26_1 = tmp10_0*w4;
                                     const register double tmp43_1 = tmp17_0*w5;  
                                     const register double tmp15_1 = A_01_2*w4;  
                                     const register double tmp53_1 = tmp19_0*w2;  
2535                                      const register double tmp27_1 = tmp3_0*w11;                                      const register double tmp27_1 = tmp3_0*w11;
2536                                        const register double tmp28_1 = tmp9_0*w7;
2537                                        const register double tmp29_1 = tmp1_0*w14;
2538                                        const register double tmp30_1 = tmp12_0*w13;
2539                                        const register double tmp31_1 = tmp11_0*w12;
2540                                      const register double tmp32_1 = tmp13_0*w2;                                      const register double tmp32_1 = tmp13_0*w2;
2541                                      const register double tmp10_1 = tmp5_0*w9;                                      const register double tmp33_1 = tmp14_0*w5;
2542                                        const register double tmp34_1 = tmp15_0*w8;
2543                                        const register double tmp35_1 = A_01_0*w8;
2544                                        const register double tmp36_1 = tmp16_0*w1;
2545                                      const register double tmp37_1 = A_10_1*w4;                                      const register double tmp37_1 = A_10_1*w4;
2546                                      const register double tmp38_1 = tmp5_0*w15;                                      const register double tmp38_1 = tmp5_0*w15;
2547                                      const register double tmp17_1 = A_01_1*w7;                                      const register double tmp39_1 = A_10_2*w7;
2548                                      const register double tmp12_1 = tmp7_0*w4;                                      const register double tmp40_1 = tmp11_0*w17;
2549                                      const register double tmp22_1 = tmp10_0*w7;                                      const register double tmp41_1 = A_01_3*w2;
                                     const register double tmp57_1 = tmp18_0*w2;  
                                     const register double tmp28_1 = tmp9_0*w7;  
                                     const register double tmp29_1 = tmp1_0*w14;  
                                     const register double tmp51_1 = tmp11_0*w16;  
2550                                      const register double tmp42_1 = tmp12_0*w16;                                      const register double tmp42_1 = tmp12_0*w16;
2551                                      const register double tmp49_1 = tmp12_0*w17;                                      const register double tmp43_1 = tmp17_0*w5;
2552                                      const register double tmp21_1 = tmp1_0*w11;                                      const register double tmp44_1 = tmp7_0*w7;
                                     const register double tmp1_1 = tmp0_0*w1;  
2553                                      const register double tmp45_1 = tmp6_0*w4;                                      const register double tmp45_1 = tmp6_0*w4;
2554                                      const register double tmp7_1 = A_10_0*w2;                                      const register double tmp46_1 = A_01_3*w8;
2555                                      const register double tmp6_1 = tmp3_0*w6;                                      const register double tmp47_1 = A_10_2*w4;
2556                                      const register double tmp13_1 = tmp8_0*w1;                                      const register double tmp48_1 = A_10_1*w7;
2557                                      const register double tmp36_1 = tmp16_0*w1;                                      const register double tmp49_1 = tmp12_0*w17;
2558                                      const register double tmp41_1 = A_01_3*w2;                                      const register double tmp50_1 = A_01_0*w2;
2559                                      const register double tmp30_1 = tmp12_0*w13;                                      const register double tmp51_1 = tmp11_0*w16;
2560                                      const register double tmp4_1 = A_01_2*w7;                                      const register double tmp52_1 = tmp18_0*w8;
2561                                      const register double tmp0_1 = A_10_3*w8;                                      const register double tmp53_1 = tmp19_0*w2;
2562                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;                                      const register double tmp54_1 = tmp13_0*w8;
2563                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;                                      const register double tmp55_1 = tmp15_0*w2;
2564                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;                                      const register double tmp56_1 = tmp19_0*w8;
2565                                        const register double tmp57_1 = tmp18_0*w2;
2566                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
2567                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;
2568                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;
2569                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;
2570                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
2571                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;
2572                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;
2573                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;
2574                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;
2575                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
2576                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;
2577                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;
2578                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;
2579                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;
2580                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;
2581                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
2582                                  }                                  }
2583                              }                              }
2584                          } else { /* constant data */                          } else { // constant data
2585                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2586                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2587                                      const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];                                      const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
# Line 2639  void Rectangle::assemblePDESystem(Paso_S Line 2590  void Rectangle::assemblePDESystem(Paso_S
2590                                      const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];                                      const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
2591                                      const register double tmp0_0 = A_01 + A_10;                                      const register double tmp0_0 = A_01 + A_10;
2592                                      const register double tmp0_1 = A_00*w18;                                      const register double tmp0_1 = A_00*w18;
2593                                      const register double tmp10_1 = A_01*w20;                                      const register double tmp1_1 = A_01*w19;
2594                                      const register double tmp12_1 = A_00*w26;                                      const register double tmp2_1 = A_10*w20;
2595                                        const register double tmp3_1 = A_11*w21;
2596                                      const register double tmp4_1 = A_00*w22;                                      const register double tmp4_1 = A_00*w22;
2597                                        const register double tmp5_1 = tmp0_0*w19;
2598                                        const register double tmp6_1 = A_11*w23;
2599                                        const register double tmp7_1 = A_11*w25;
2600                                      const register double tmp8_1 = A_00*w24;                                      const register double tmp8_1 = A_00*w24;
                                     const register double tmp13_1 = A_10*w19;  
2601                                      const register double tmp9_1 = tmp0_0*w20;                                      const register double tmp9_1 = tmp0_0*w20;
2602                                      const register double tmp3_1 = A_11*w21;                                      const register double tmp10_1 = A_01*w20;
2603                                      const register double tmp11_1 = A_11*w27;                                      const register double tmp11_1 = A_11*w27;
2604                                      const register double tmp1_1 = A_01*w19;                                      const register double tmp12_1 = A_00*w26;
2605                                      const register double tmp6_1 = A_11*w23;                                      const register double tmp13_1 = A_10*w19;
                                     const register double tmp7_1 = A_11*w25;  
                                     const register double tmp2_1 = A_10*w20;  
                                     const register double tmp5_1 = tmp0_0*w19;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
2606                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2607                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2608                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
2609                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2610                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2611                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2612                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
2613                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
2614                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
2615                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
2616                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
2617                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2618                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;
2619                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;
2620                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;
2621                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;
2622                                  }                                  }
2623                              }                              }
2624                          }                          }
# Line 2689  void Rectangle::assemblePDESystem(Paso_S Line 2640  void Rectangle::assemblePDESystem(Paso_S
2640                                      const register double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];                                      const register double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
2641                                      const register double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];                                      const register double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
2642                                      const register double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];                                      const register double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
2643                                      const register double tmp3_0 = B_0_0 + B_0_2;                                      const register double tmp0_0 = B_1_0 + B_1_1;
2644                                      const register double tmp1_0 = B_1_2 + B_1_3;                                      const register double tmp1_0 = B_1_2 + B_1_3;
2645                                      const register double tmp2_0 = B_0_1 + B_0_3;                                      const register double tmp2_0 = B_0_1 + B_0_3;
2646                                      const register double tmp0_0 = B_1_0 + B_1_1;                                      const register double tmp3_0 = B_0_0 + B_0_2;
2647                                      const register double tmp63_1 = B_1_1*w42;                                      const register double tmp63_1 = B_1_1*w42;
2648                                      const register double tmp79_1 = B_1_1*w40;                                      const register double tmp79_1 = B_1_1*w40;
2649                                      const register double tmp37_1 = tmp3_0*w35;                                      const register double tmp37_1 = tmp3_0*w35;
# Line 2791  void Rectangle::assemblePDESystem(Paso_S Line 2742  void Rectangle::assemblePDESystem(Paso_S
2742                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
2743                                  }                                  }
2744                              }                              }
2745                          } else { /* constant data */                          } else { // constant data
2746                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2747                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2748                                      const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];                                      const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
# Line 2943  void Rectangle::assemblePDESystem(Paso_S Line 2894  void Rectangle::assemblePDESystem(Paso_S
2894                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;
2895                                  }                                  }
2896                              }                              }
2897                          } else { /* constant data */                          } else { // constant data
2898                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2899                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2900                                      const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];                                      const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
# Line 3033  void Rectangle::assemblePDESystem(Paso_S Line 2984  void Rectangle::assemblePDESystem(Paso_S
2984                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp14_1 + tmp17_1 + tmp18_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp14_1 + tmp17_1 + tmp18_1;
2985                                  }                                  }
2986                               }                               }
2987                          } else { /* constant data */                          } else { // constant data
2988                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2989                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2990                                      const register double D_0 = D_p[INDEX2(k, m, numEq)];                                      const register double D_0 = D_p[INDEX2(k, m, numEq)];
# Line 3101  void Rectangle::assemblePDESystem(Paso_S Line 3052  void Rectangle::assemblePDESystem(Paso_S
3052                                  EM_F[INDEX2(k,2,numEq)]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;
3053                                  EM_F[INDEX2(k,3,numEq)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;                                  EM_F[INDEX2(k,3,numEq)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
3054                              }                              }
3055                          } else { /* constant data */                          } else { // constant data
3056                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3057                                  const register double X_0 = X_p[INDEX2(k, 0, numEq)];                                  const register double X_0 = X_p[INDEX2(k, 0, numEq)];
3058                                  const register double X_1 = X_p[INDEX2(k, 1, numEq)];                                  const register double X_1 = X_p[INDEX2(k, 1, numEq)];
                                 const register double tmp3_1 = X_1*w71;  
                                 const register double tmp2_1 = X_0*w70;  
                                 const register double tmp1_1 = X_0*w68;  
3059                                  const register double tmp0_1 = X_1*w69;                                  const register double tmp0_1 = X_1*w69;
3060                                    const register double tmp1_1 = X_0*w68;
3061                                    const register double tmp2_1 = X_0*w70;
3062                                    const register double tmp3_1 = X_1*w71;
3063                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3064                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;
3065                                  EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;
# Line 3130  void Rectangle::assemblePDESystem(Paso_S Line 3081  void Rectangle::assemblePDESystem(Paso_S
3081                                  const register double Y_3 = Y_p[INDEX2(k, 3, numEq)];                                  const register double Y_3 = Y_p[INDEX2(k, 3, numEq)];
3082                                  const register double tmp0_0 = Y_1 + Y_2;                                  const register double tmp0_0 = Y_1 + Y_2;
3083                                  const register double tmp1_0 = Y_0 + Y_3;                                  const register double tmp1_0 = Y_0 + Y_3;
3084                                  const register double tmp9_1 = Y_0*w74;                                  const register double tmp0_1 = Y_0*w72;
3085                                    const register double tmp1_1 = tmp0_0*w73;
3086                                    const register double tmp2_1 = Y_3*w74;
3087                                    const register double tmp3_1 = Y_1*w72;
3088                                  const register double tmp4_1 = tmp1_0*w73;                                  const register double tmp4_1 = tmp1_0*w73;
3089                                  const register double tmp5_1 = Y_2*w74;                                  const register double tmp5_1 = Y_2*w74;
                                 const register double tmp7_1 = Y_1*w74;  
3090                                  const register double tmp6_1 = Y_2*w72;                                  const register double tmp6_1 = Y_2*w72;
3091                                  const register double tmp2_1 = Y_3*w74;                                  const register double tmp7_1 = Y_1*w74;
3092                                  const register double tmp8_1 = Y_3*w72;                                  const register double tmp8_1 = Y_3*w72;
3093                                  const register double tmp3_1 = Y_1*w72;                                  const register double tmp9_1 = Y_0*w74;
                                 const register double tmp0_1 = Y_0*w72;  
                                 const register double tmp1_1 = tmp0_0*w73;  
3094                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
3095                                  EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
3096                                  EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;
3097                                  EM_F[INDEX2(k,3,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;                                  EM_F[INDEX2(k,3,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;
3098                              }                              }
3099                          } else { /* constant data */                          } else { // constant data
3100                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3101                                  const register double Y_0 = Y_p[k];                                  const register double tmp0_1 = Y_p[k]*w75;
                                 const register double tmp0_1 = Y_0*w75;  
3102                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3103                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3104                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
# Line 3156  void Rectangle::assemblePDESystem(Paso_S Line 3106  void Rectangle::assemblePDESystem(Paso_S
3106                              }                              }
3107                          }                          }
3108                      }                      }
                     /* GENERATOR SNIP_PDE_SYSTEM BOTTOM */  
3109    
3110                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
3111                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
# Line 3166  void Rectangle::assemblePDESystem(Paso_S Line 3115  void Rectangle::assemblePDESystem(Paso_S
3115                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);
3116                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
3117                      if (add_EM_F) {                      if (add_EM_F) {
                         //cout << "-- AddtoRHS -- " << endl;  
3118                          double *F_p=rhs.getSampleDataRW(0);                          double *F_p=rhs.getSampleDataRW(0);
3119                          for (index_t i=0; i<rowIndex.size(); i++) {                          for (index_t i=0; i<rowIndex.size(); i++) {
3120                              if (rowIndex[i]<getNumDOF()) {                              if (rowIndex[i]<getNumDOF()) {
3121                                  for (index_t eq=0; eq<numEq; eq++) {                                  for (index_t eq=0; eq<numEq; eq++) {
3122                                      F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];                                      F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];
                                     //cout << "F[" << INDEX2(eq,rowIndex[i],numEq) << "]=" << F_p[INDEX2(eq,rowIndex[i],numEq)] << endl;  
3123                                  }                                  }
3124                              }                              }
3125                          }                          }
                         //cout << "---"<<endl;  
3126                      }                      }
3127                      if (add_EM_S) {                      if (add_EM_S) {
                         //cout << "-- AddtoSystem -- " << endl;  
3128                          addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);                          addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);
3129                      }                      }
3130    
# Line 3206  void Rectangle::assemblePDESystemReduced Line 3151  void Rectangle::assemblePDESystemReduced
3151          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
3152      }      }
3153    
3154      /*** GENERATOR SNIP_PDE_SYSTEM_REDUCED_PRE TOP */      const double w0 = -.25*h1/h0;
3155      const double w0 = -0.25*h1/h0;      const double w1 = .25;
3156      const double w1 = 0.25;      const double w2 = -.25;
3157      const double w2 = -0.25;      const double w3 = .25*h0/h1;
3158      const double w3 = 0.25*h0/h1;      const double w4 = -.25*h0/h1;
3159      const double w4 = -0.25*h0/h1;      const double w5 = .25*h1/h0;
3160      const double w5 = 0.25*h1/h0;      const double w6 = -.125*h1;
3161      const double w6 = -0.125*h1;      const double w7 = -.125*h0;
3162      const double w7 = -0.125*h0;      const double w8 = .125*h1;
3163      const double w8 = 0.125*h1;      const double w9 = .125*h0;
3164      const double w9 = 0.125*h0;      const double w10 = .0625*h0*h1;
3165      const double w10 = 0.0625*h0*h1;      const double w11 = -.5*h1;
3166      const double w11 = -0.5*h1;      const double w12 = -.5*h0;
3167      const double w12 = -0.5*h0;      const double w13 = .5*h1;
3168      const double w13 = 0.5*h1;      const double w14 = .5*h0;
3169      const double w14 = 0.5*h0;      const double w15 = .25*h0*h1;
     const double w15 = 0.25*h0*h1;  
     /* GENERATOR SNIP_PDE_SYSTEM_REDUCED_PRE BOTTOM */  
3170    
3171      rhs.requireWrite();      rhs.requireWrite();
3172  #pragma omp parallel  #pragma omp parallel
# Line 3237  void Rectangle::assemblePDESystemReduced Line 3180  void Rectangle::assemblePDESystemReduced
3180                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
3181                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
3182                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
                     /*** GENERATOR SNIP_PDE_SYSTEM_REDUCED TOP */  
3183                      ///////////////                      ///////////////
3184                      // process A //                      // process A //
3185                      ///////////////                      ///////////////
# Line 3251  void Rectangle::assemblePDESystemReduced Line 3193  void Rectangle::assemblePDESystemReduced
3193                                  const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];                                  const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
3194                                  const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];                                  const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
3195                                  const register double tmp0_0 = A_01 + A_10;                                  const register double tmp0_0 = A_01 + A_10;
3196                                    const register double tmp0_1 = A_11*w3;
3197                                    const register double tmp1_1 = A_00*w0;
3198                                  const register double tmp2_1 = A_01*w1;                                  const register double tmp2_1 = A_01*w1;
                                 const register double tmp7_1 = tmp0_0*w2;  
3199                                  const register double tmp3_1 = A_10*w2;                                  const register double tmp3_1 = A_10*w2;
                                 const register double tmp6_1 = A_00*w5;  
                                 const register double tmp1_1 = A_00*w0;  
3200                                  const register double tmp4_1 = tmp0_0*w1;                                  const register double tmp4_1 = tmp0_0*w1;
                                 const register double tmp9_1 = A_01*w2;  
3201                                  const register double tmp5_1 = A_11*w4;                                  const register double tmp5_1 = A_11*w4;
3202                                    const register double tmp6_1 = A_00*w5;
3203                                    const register double tmp7_1 = tmp0_0*w2;
3204                                  const register double tmp8_1 = A_10*w1;                                  const register double tmp8_1 = A_10*w1;
3205                                  const register double tmp0_1 = A_11*w3;                                  const register double tmp9_1 = A_01*w2;
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
3206                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
3207                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
3208                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
3209                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
3210                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3211                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
3212                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
3213                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
3214                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
3215                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
3216                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
3217                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3218                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
3219                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
3220                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
3221                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
3222                              }                              }
3223                          }                          }
3224                      }                      }
# Line 3290  void Rectangle::assemblePDESystemReduced Line 3232  void Rectangle::assemblePDESystemReduced
3232                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3233                                  const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];                                  const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
3234                                  const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];                                  const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
                                 const register double tmp2_1 = B_0*w8;  
3235                                  const register double tmp0_1 = B_0*w6;                                  const register double tmp0_1 = B_0*w6;
                                 const register double tmp3_1 = B_1*w9;  
3236                                  const register double tmp1_1 = B_1*w7;                                  const register double tmp1_1 = B_1*w7;
3237                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  const register double tmp2_1 = B_0*w8;
3238                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                  const register double tmp3_1 = B_1*w9;
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
3239                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3240                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3241                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3242                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3243                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3244                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3245                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3246                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3247                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3248                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp3_1;  
3249                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3250                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3251                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3252                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3253                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3254                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3255                              }                              }
3256                          }                          }
3257                      }                      }
# Line 3323  void Rectangle::assemblePDESystemReduced Line 3265  void Rectangle::assemblePDESystemReduced
3265                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3266                                  const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];                                  const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
3267                                  const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];                                  const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
                                 const register double tmp1_1 = C_1*w7;  
3268                                  const register double tmp0_1 = C_0*w8;                                  const register double tmp0_1 = C_0*w8;
3269                                  const register double tmp3_1 = C_0*w6;                                  const register double tmp1_1 = C_1*w7;
3270                                  const register double tmp2_1 = C_1*w9;                                  const register double tmp2_1 = C_1*w9;
3271                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  const register double tmp3_1 = C_0*w6;
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
3272                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3273                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3274                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3275                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3276                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3277                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3278                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3279                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3280                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3281                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
3282                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3283                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3284                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3285                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3286                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3287                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3288                              }                              }
3289                          }                          }
3290                      }                      }
# Line 3354  void Rectangle::assemblePDESystemReduced Line 3296  void Rectangle::assemblePDESystemReduced
3296                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
3297                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3298                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3299                                  const register double D_0 = D_p[INDEX2(k, m, numEq)];                                  const register double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;
                                 const register double tmp0_1 = D_0*w10;  
                                 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;  
3300                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
3301                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
3302                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;
3303                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;
3304                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
3305                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
3306                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1;
3307                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;
3308                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;
3309                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;  
3310                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
3311                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
3312                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1;
3313                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;
3314                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
3315                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
3316                              }                              }
3317                          }                          }
3318                      }                      }
# Line 3385  void Rectangle::assemblePDESystemReduced Line 3326  void Rectangle::assemblePDESystemReduced
3326                              const register double X_0 = X_p[INDEX2(k, 0, numEq)];                              const register double X_0 = X_p[INDEX2(k, 0, numEq)];
3327                              const register double X_1 = X_p[INDEX2(k, 1, numEq)];                              const register double X_1 = X_p[INDEX2(k, 1, numEq)];
3328                              const register double tmp0_1 = X_0*w11;                              const register double tmp0_1 = X_0*w11;
                             const register double tmp2_1 = X_0*w13;  
3329                              const register double tmp1_1 = X_1*w12;                              const register double tmp1_1 = X_1*w12;
3330                                const register double tmp2_1 = X_0*w13;
3331                              const register double tmp3_1 = X_1*w14;                              const register double tmp3_1 = X_1*w14;
3332                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3333                              EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;
# Line 3401  void Rectangle::assemblePDESystemReduced Line 3342  void Rectangle::assemblePDESystemReduced
3342                          add_EM_F=true;                          add_EM_F=true;
3343                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3344                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3345                              const register double Y_0 = Y_p[k];                              const register double tmp0_1 = Y_p[k]*w15;
                             const register double tmp0_1 = Y_0*w15;  
3346                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3347                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3348                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
3349                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
3350                          }                          }
3351                      }                      }
                     /* GENERATOR SNIP_PDE_SYSTEM_REDUCED BOTTOM */  
3352    
3353                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
3354                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
# Line 3419  void Rectangle::assemblePDESystemReduced Line 3358  void Rectangle::assemblePDESystemReduced
3358                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);
3359                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
3360                      if (add_EM_F) {                      if (add_EM_F) {
                         //cout << "-- AddtoRHS -- " << endl;  
3361                          double *F_p=rhs.getSampleDataRW(0);                          double *F_p=rhs.getSampleDataRW(0);
3362                          for (index_t i=0; i<rowIndex.size(); i++) {                          for (index_t i=0; i<rowIndex.size(); i++) {
3363                              if (rowIndex[i]<getNumDOF()) {                              if (rowIndex[i]<getNumDOF()) {
3364                                  for (index_t eq=0; eq<numEq; eq++) {                                  for (index_t eq=0; eq<numEq; eq++) {
3365                                      F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];                                      F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];
                                     //cout << "F[" << INDEX2(eq,rowIndex[i],numEq) << "]=" << F_p[INDEX2(eq,rowIndex[i],numEq)] << endl;  
3366                                  }                                  }
3367                              }                              }
3368                          }                          }
                         //cout << "---"<<endl;  
3369                      }                      }
3370                      if (add_EM_S) {                      if (add_EM_S) {
                         //cout << "-- AddtoSystem -- " << endl;  
3371                          addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);                          addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);
3372                      }                      }
3373    

Legend:
Removed from v.3762  
changed lines
  Added in v.3764

  ViewVC Help
Powered by ViewVC 1.1.26