/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Diff of /branches/diaplayground/ripley/src/Brick.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 54  Brick::Brick(int n0, int n1, int n2, dou Line 54  Brick::Brick(int n0, int n1, int n2, dou
54          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
55    
56      // local number of elements (including overlap)      // local number of elements (including overlap)
57      m_NE0 = (m_NX>1 ? (n0+1)/m_NX : n0);      m_NE0 = m_ownNE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
58      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)
59          m_NE0++;          m_NE0++;
60      m_NE1 = (m_NY>1 ? (n1+1)/m_NY : n1);      else if (m_NX>1 && m_mpiInfo->rank%m_NX==m_NX-1)
61            m_ownNE0--;
62    
63        m_NE1 = m_ownNE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
64      if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX>0 && m_mpiInfo->rank%(m_NX*m_NY)/m_NX<m_NY-1)      if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX>0 && m_mpiInfo->rank%(m_NX*m_NY)/m_NX<m_NY-1)
65          m_NE1++;          m_NE1++;
66      m_NE2 = (m_NZ>1 ? (n2+1)/m_NZ : n2);      else if (m_NY>1 && m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
67            m_ownNE1--;
68    
69        m_NE2 = m_ownNE2 = (m_NZ>1 ? (n2+1)/m_NZ : n2);
70      if (m_mpiInfo->rank/(m_NX*m_NY)>0 && m_mpiInfo->rank/(m_NX*m_NY)<m_NZ-1)      if (m_mpiInfo->rank/(m_NX*m_NY)>0 && m_mpiInfo->rank/(m_NX*m_NY)<m_NZ-1)
71          m_NE2++;          m_NE2++;
72        else if (m_NZ>1 && m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
73            m_ownNE2--;
74    
75      // local number of nodes      // local number of nodes
76      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
# Line 332  bool Brick::ownSample(int fsType, index_ Line 340  bool Brick::ownSample(int fsType, index_
340      throw RipleyException(msg.str());      throw RipleyException(msg.str());
341  }  }
342    
343  void Brick::setToGradient(escript::Data& out, const escript::Data& cIn) const  void Brick::setToNormal(escript::Data& out) const
344    {
345        if (out.getFunctionSpace().getTypeCode() == FaceElements) {
346            out.requireWrite();
347    #pragma omp parallel
348            {
349                if (m_faceOffset[0] > -1) {
350    #pragma omp for nowait
351                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
352                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
353                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
354                            // set vector at four quadrature points
355                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
356                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
357                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
358                            *o++ = -1.; *o++ = 0.; *o = 0.;
359                        }
360                    }
361                }
362    
363                if (m_faceOffset[1] > -1) {
364    #pragma omp for nowait
365                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
366                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
367                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
368                            // set vector at four quadrature points
369                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
370                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
371                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
372                            *o++ = 1.; *o++ = 0.; *o = 0.;
373                        }
374                    }
375                }
376    
377                if (m_faceOffset[2] > -1) {
378    #pragma omp for nowait
379                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
380                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
381                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
382                            // set vector at four quadrature points
383                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
384                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
385                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
386                            *o++ = 0.; *o++ = -1.; *o = 0.;
387                        }
388                    }
389                }
390    
391                if (m_faceOffset[3] > -1) {
392    #pragma omp for nowait
393                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
394                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
395                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
396                            // set vector at four quadrature points
397                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
398                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
399                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
400                            *o++ = 0.; *o++ = 1.; *o = 0.;
401                        }
402                    }
403                }
404    
405                if (m_faceOffset[4] > -1) {
406    #pragma omp for nowait
407                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
408                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
409                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
410                            // set vector at four quadrature points
411                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
412                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
413                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
414                            *o++ = 0.; *o++ = 0.; *o = -1.;
415                        }
416                    }
417                }
418    
419                if (m_faceOffset[5] > -1) {
420    #pragma omp for nowait
421                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
422                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
423                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
424                            // set vector at four quadrature points
425                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
426                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
427                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
428                            *o++ = 0.; *o++ = 0.; *o = 1.;
429                        }
430                    }
431                }
432            } // end of parallel section
433        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
434            out.requireWrite();
435    #pragma omp parallel
436            {
437                if (m_faceOffset[0] > -1) {
438    #pragma omp for nowait
439                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
440                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
441                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
442                            *o++ = -1.;
443                            *o++ = 0.;
444                            *o = 0.;
445                        }
446                    }
447                }
448    
449                if (m_faceOffset[1] > -1) {
450    #pragma omp for nowait
451                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
452                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
453                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
454                            *o++ = 1.;
455                            *o++ = 0.;
456                            *o = 0.;
457                        }
458                    }
459                }
460    
461                if (m_faceOffset[2] > -1) {
462    #pragma omp for nowait
463                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
464                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
465                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
466                            *o++ = 0.;
467                            *o++ = -1.;
468                            *o = 0.;
469                        }
470                    }
471                }
472    
473                if (m_faceOffset[3] > -1) {
474    #pragma omp for nowait
475                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
476                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
477                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
478                            *o++ = 0.;
479                            *o++ = 1.;
480                            *o = 0.;
481                        }
482                    }
483                }
484    
485                if (m_faceOffset[4] > -1) {
486    #pragma omp for nowait
487                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
488                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
489                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
490                            *o++ = 0.;
491                            *o++ = 0.;
492                            *o = -1.;
493                        }
494                    }
495                }
496    
497                if (m_faceOffset[5] > -1) {
498    #pragma omp for nowait
499                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
500                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
501                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
502                            *o++ = 0.;
503                            *o++ = 0.;
504                            *o = 1.;
505                        }
506                    }
507                }
508            } // end of parallel section
509    
510        } else {
511            stringstream msg;
512            msg << "setToNormal() not implemented for "
513                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
514            throw RipleyException(msg.str());
515        }
516    }
517    
518    void Brick::setToSize(escript::Data& out) const
519    {
520        if (out.getFunctionSpace().getTypeCode() == Elements
521                || out.getFunctionSpace().getTypeCode() == ReducedElements) {
522            out.requireWrite();
523            const dim_t numQuad=out.getNumDataPointsPerSample();
524            const double xSize=getFirstCoordAndSpacing(0).second;
525            const double ySize=getFirstCoordAndSpacing(1).second;
526            const double zSize=getFirstCoordAndSpacing(2).second;
527            const double size=min(min(xSize,ySize),zSize);
528    #pragma omp parallel for
529            for (index_t k = 0; k < getNumElements(); ++k) {
530                double* o = out.getSampleDataRW(k);
531                fill(o, o+numQuad, size);
532            }
533        } else if (out.getFunctionSpace().getTypeCode() == FaceElements
534                || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
535            out.requireWrite();
536            const dim_t numQuad=out.getNumDataPointsPerSample();
537            const double xSize=getFirstCoordAndSpacing(0).second;
538            const double ySize=getFirstCoordAndSpacing(1).second;
539            const double zSize=getFirstCoordAndSpacing(2).second;
540    #pragma omp parallel
541            {
542                if (m_faceOffset[0] > -1) {
543                    const double size=min(ySize,zSize);
544    #pragma omp for nowait
545                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
546                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
547                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
548                            fill(o, o+numQuad, size);
549                        }
550                    }
551                }
552    
553                if (m_faceOffset[1] > -1) {
554                    const double size=min(ySize,zSize);
555    #pragma omp for nowait
556                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
557                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
558                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
559                            fill(o, o+numQuad, size);
560                        }
561                    }
562                }
563    
564                if (m_faceOffset[2] > -1) {
565                    const double size=min(xSize,zSize);
566    #pragma omp for nowait
567                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
568                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
569                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
570                            fill(o, o+numQuad, size);
571                        }
572                    }
573                }
574    
575                if (m_faceOffset[3] > -1) {
576                    const double size=min(xSize,zSize);
577    #pragma omp for nowait
578                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
579                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
580                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
581                            fill(o, o+numQuad, size);
582                        }
583                    }
584                }
585    
586                if (m_faceOffset[4] > -1) {
587                    const double size=min(xSize,ySize);
588    #pragma omp for nowait
589                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
590                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
591                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
592                            fill(o, o+numQuad, size);
593                        }
594                    }
595                }
596    
597                if (m_faceOffset[5] > -1) {
598                    const double size=min(xSize,ySize);
599    #pragma omp for nowait
600                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
601                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
602                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
603                            fill(o, o+numQuad, size);
604                        }
605                    }
606                }
607            } // end of parallel section
608    
609        } else {
610            stringstream msg;
611            msg << "setToSize() not implemented for "
612                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
613            throw RipleyException(msg.str());
614        }
615    }
616    
617    Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
618                                                bool reducedColOrder) const
619    {
620        if (reducedRowOrder || reducedColOrder)
621            throw RipleyException("getPattern() not implemented for reduced order");
622    
623        return m_pattern;
624    }
625    
626    void Brick::Print_Mesh_Info(const bool full) const
627    {
628        RipleyDomain::Print_Mesh_Info(full);
629        if (full) {
630            cout << "     Id  Coordinates" << endl;
631            cout.precision(15);
632            cout.setf(ios::scientific, ios::floatfield);
633            pair<double,double> xdx = getFirstCoordAndSpacing(0);
634            pair<double,double> ydy = getFirstCoordAndSpacing(1);
635            pair<double,double> zdz = getFirstCoordAndSpacing(2);
636            for (index_t i=0; i < getNumNodes(); i++) {
637                cout << "  " << setw(5) << m_nodeId[i]
638                    << "  " << xdx.first+(i%m_N0)*xdx.second
639                    << "  " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second
640                    << "  " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;
641            }
642        }
643    }
644    
645    IndexVector Brick::getNumNodesPerDim() const
646    {
647        IndexVector ret;
648        ret.push_back(m_N0);
649        ret.push_back(m_N1);
650        ret.push_back(m_N2);
651        return ret;
652    }
653    
654    IndexVector Brick::getNumElementsPerDim() const
655    {
656        IndexVector ret;
657        ret.push_back(m_NE0);
658        ret.push_back(m_NE1);
659        ret.push_back(m_NE2);
660        return ret;
661    }
662    
663    IndexVector Brick::getNumFacesPerBoundary() const
664    {
665        IndexVector ret(6, 0);
666        //left
667        if (m_offset0==0)
668            ret[0]=m_NE1*m_NE2;
669        //right
670        if (m_mpiInfo->rank%m_NX==m_NX-1)
671            ret[1]=m_NE1*m_NE2;
672        //bottom
673        if (m_offset1==0)
674            ret[2]=m_NE0*m_NE2;
675        //top
676        if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
677            ret[3]=m_NE0*m_NE2;
678        //front
679        if (m_offset2==0)
680            ret[4]=m_NE0*m_NE1;
681        //back
682        if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
683            ret[5]=m_NE0*m_NE1;
684        return ret;
685    }
686    
687    pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const
688    {
689        if (dim==0)
690            return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
691        else if (dim==1)
692            return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
693        else if (dim==2)
694            return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);
695    
696        throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
697    }
698    
699    //protected
700    dim_t Brick::getNumDOF() const
701    {
702        return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY*(m_gNE2+1)/m_NZ;
703    }
704    
705    //protected
706    dim_t Brick::getNumFaceElements() const
707    {
708        const IndexVector faces = getNumFacesPerBoundary();
709        dim_t n=0;
710        for (size_t i=0; i<faces.size(); i++)
711            n+=faces[i];
712        return n;
713    }
714    
715    //protected
716    void Brick::assembleCoordinates(escript::Data& arg) const
717    {
718        escriptDataC x = arg.getDataC();
719        int numDim = m_numDim;
720        if (!isDataPointShapeEqual(&x, 1, &numDim))
721            throw RipleyException("setToX: Invalid Data object shape");
722        if (!numSamplesEqual(&x, 1, getNumNodes()))
723            throw RipleyException("setToX: Illegal number of samples in Data object");
724    
725        pair<double,double> xdx = getFirstCoordAndSpacing(0);
726        pair<double,double> ydy = getFirstCoordAndSpacing(1);
727        pair<double,double> zdz = getFirstCoordAndSpacing(2);
728        arg.requireWrite();
729    #pragma omp parallel for
730        for (dim_t i2 = 0; i2 < m_N2; i2++) {
731            for (dim_t i1 = 0; i1 < m_N1; i1++) {
732                for (dim_t i0 = 0; i0 < m_N0; i0++) {
733                    double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);
734                    point[0] = xdx.first+i0*xdx.second;
735                    point[1] = ydy.first+i1*ydy.second;
736                    point[2] = zdz.first+i2*zdz.second;
737                }
738            }
739        }
740    }
741    
742    //protected
743    void Brick::assembleGradient(escript::Data& out, escript::Data& in) const
744  {  {
     escript::Data& in = *const_cast<escript::Data*>(&cIn);  
745      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
746      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
747      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
# Line 349  void Brick::setToGradient(escript::Data& Line 756  void Brick::setToGradient(escript::Data&
756    
757      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
758          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */  
759  #pragma omp parallel for  #pragma omp parallel for
760          for (index_t k2=0; k2 < m_NE2; ++k2) {          for (index_t k2=0; k2 < m_NE2; ++k2) {
761              for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 400  void Brick::setToGradient(escript::Data& Line 806  void Brick::setToGradient(escript::Data&
806                          o[INDEX3(i,0,7,numComp,3)] = V3;                          o[INDEX3(i,0,7,numComp,3)] = V3;
807                          o[INDEX3(i,1,7,numComp,3)] = V7;                          o[INDEX3(i,1,7,numComp,3)] = V7;
808                          o[INDEX3(i,2,7,numComp,3)] = V11;                          o[INDEX3(i,2,7,numComp,3)] = V11;
809                      } /* end of component loop i */                      } // end of component loop i
810                  } /* end of k0 loop */                  } // end of k0 loop
811              } /* end of k1 loop */              } // end of k1 loop
812          } /* end of k2 loop */          } // end of k2 loop
         /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */  
813      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
814          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */  
815  #pragma omp parallel for  #pragma omp parallel for
816          for (index_t k2=0; k2 < m_NE2; ++k2) {          for (index_t k2=0; k2 < m_NE2; ++k2) {
817              for (index_t k1=0; k1 < m_NE1; ++k1) {              for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 425  void Brick::setToGradient(escript::Data& Line 829  void Brick::setToGradient(escript::Data&
829                          o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;                          o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;
830                          o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;                          o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;
831                          o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / h2;                          o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / h2;
832                      } /* end of component loop i */                      } // end of component loop i
833                  } /* end of k0 loop */                  } // end of k0 loop
834              } /* end of k1 loop */              } // end of k1 loop
835          } /* end of k2 loop */          } // end of k2 loop
         /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */  
836      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
837          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_GRAD_FACES TOP */  
838  #pragma omp parallel  #pragma omp parallel
839          {          {
840              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
# Line 465  void Brick::setToGradient(escript::Data& Line 867  void Brick::setToGradient(escript::Data&
867                              o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;                              o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
868                              o[INDEX3(i,1,3,numComp,3)] = V1;                              o[INDEX3(i,1,3,numComp,3)] = V1;
869                              o[INDEX3(i,2,3,numComp,3)] = V3;                              o[INDEX3(i,2,3,numComp,3)] = V3;
870                          } /* end of component loop i */                          } // end of component loop i
871                      } /* end of k1 loop */                      } // end of k1 loop
872                  } /* end of k2 loop */                  } // end of k2 loop
873              } /* end of face 0 */              } // end of face 0
874              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
875  #pragma omp for nowait  #pragma omp for nowait
876                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 499  void Brick::setToGradient(escript::Data& Line 901  void Brick::setToGradient(escript::Data&
901                              o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;                              o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
902                              o[INDEX3(i,1,3,numComp,3)] = V1;                              o[INDEX3(i,1,3,numComp,3)] = V1;
903                              o[INDEX3(i,2,3,numComp,3)] = V3;                              o[INDEX3(i,2,3,numComp,3)] = V3;
904                          } /* end of component loop i */                          } // end of component loop i
905                      } /* end of k1 loop */                      } // end of k1 loop
906                  } /* end of k2 loop */                  } // end of k2 loop
907              } /* end of face 1 */              } // end of face 1
908              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
909  #pragma omp for nowait  #pragma omp for nowait
910                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 532  void Brick::setToGradient(escript::Data& Line 934  void Brick::setToGradient(escript::Data&
934                              o[INDEX3(i,0,3,numComp,3)] = V0;                              o[INDEX3(i,0,3,numComp,3)] = V0;
935                              o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;                              o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
936                              o[INDEX3(i,2,3,numComp,3)] = V2;                              o[INDEX3(i,2,3,numComp,3)] = V2;
937                          } /* end of component loop i */                          } // end of component loop i
938                      } /* end of k0 loop */                      } // end of k0 loop
939                  } /* end of k2 loop */                  } // end of k2 loop
940              } /* end of face 2 */              } // end of face 2
941              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
942  #pragma omp for nowait  #pragma omp for nowait
943                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 566  void Brick::setToGradient(escript::Data& Line 968  void Brick::setToGradient(escript::Data&
968                              o[INDEX3(i,0,3,numComp,3)] = V1;                              o[INDEX3(i,0,3,numComp,3)] = V1;
969                              o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;                              o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
970                              o[INDEX3(i,2,3,numComp,3)] = V3;                              o[INDEX3(i,2,3,numComp,3)] = V3;
971                          } /* end of component loop i */                          } // end of component loop i
972                      } /* end of k0 loop */                      } // end of k0 loop
973                  } /* end of k2 loop */                  } // end of k2 loop
974              } /* end of face 3 */              } // end of face 3
975              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
976  #pragma omp for nowait  #pragma omp for nowait
977                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 600  void Brick::setToGradient(escript::Data& Line 1002  void Brick::setToGradient(escript::Data&
1002                              o[INDEX3(i,0,3,numComp,3)] = V1;                              o[INDEX3(i,0,3,numComp,3)] = V1;
1003                              o[INDEX3(i,1,3,numComp,3)] = V3;                              o[INDEX3(i,1,3,numComp,3)] = V3;
1004                              o[INDEX3(i,2,3,numComp,3)] = ((f_111[i]-f_110[i])*C5 + (f_001[i]-f_000[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;                              o[INDEX3(i,2,3,numComp,3)] = ((f_111[i]-f_110[i])*C5 + (f_001[i]-f_000[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
1005                          } /* end of component loop i */                          } // end of component loop i
1006                      } /* end of k0 loop */                      } // end of k0 loop
1007                  } /* end of k1 loop */                  } // end of k1 loop
1008              } /* end of face 4 */              } // end of face 4
1009              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1010  #pragma omp for nowait  #pragma omp for nowait
1011                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 634  void Brick::setToGradient(escript::Data& Line 1036  void Brick::setToGradient(escript::Data&
1036                              o[INDEX3(i,0,3,numComp,3)] = V1;                              o[INDEX3(i,0,3,numComp,3)] = V1;
1037                              o[INDEX3(i,1,3,numComp,3)] = V3;                              o[INDEX3(i,1,3,numComp,3)] = V3;
1038                              o[INDEX3(i,2,3,numComp,3)] = ((f_001[i]-f_000[i])*C0 + (f_111[i]-f_110[i])*C5 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;                              o[INDEX3(i,2,3,numComp,3)] = ((f_001[i]-f_000[i])*C0 + (f_111[i]-f_110[i])*C5 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
1039                          } /* end of component loop i */                          } // end of component loop i
1040                      } /* end of k0 loop */                      } // end of k0 loop
1041                  } /* end of k1 loop */                  } // end of k1 loop
1042              } /* end of face 5 */              } // end of face 5
1043          } // end of parallel section          } // end of parallel section
         /* GENERATOR SNIP_GRAD_FACES BOTTOM */  
1044      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1045          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */  
1046  #pragma omp parallel  #pragma omp parallel
1047          {          {
1048              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
# Line 662  void Brick::setToGradient(escript::Data& Line 1062  void Brick::setToGradient(escript::Data&
1062                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;
1063                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / h1;
1064                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]-f_000[i]-f_010[i])*C4 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]-f_000[i]-f_010[i])*C4 / h2;
1065                          } /* end of component loop i */                          } // end of component loop i
1066                      } /* end of k1 loop */                      } // end of k1 loop
1067                  } /* end of k2 loop */                  } // end of k2 loop
1068              } /* end of face 0 */              } // end of face 0
1069              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1070  #pragma omp for nowait  #pragma omp for nowait
1071                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 683  void Brick::setToGradient(escript::Data& Line 1083  void Brick::setToGradient(escript::Data&
1083                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;
1084                              o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / h1;
1085                              o[INDEX3(i,2,0,numComp,3)] = (f_101[i]+f_111[i]-f_100[i]-f_110[i])*C4 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_101[i]+f_111[i]-f_100[i]-f_110[i])*C4 / h2;
1086                          } /* end of component loop i */                          } // end of component loop i
1087                      } /* end of k1 loop */                      } // end of k1 loop
1088                  } /* end of k2 loop */                  } // end of k2 loop
1089              } /* end of face 1 */              } // end of face 1
1090              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1091  #pragma omp for nowait  #pragma omp for nowait
1092                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 704  void Brick::setToGradient(escript::Data& Line 1104  void Brick::setToGradient(escript::Data&
1104                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / h0;
1105                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;
1106                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_101[i]-f_000[i]-f_100[i])*C4 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_101[i]-f_000[i]-f_100[i])*C4 / h2;
1107                          } /* end of component loop i */                          } // end of component loop i
1108                      } /* end of k0 loop */                      } // end of k0 loop
1109                  } /* end of k2 loop */                  } // end of k2 loop
1110              } /* end of face 2 */              } // end of face 2
1111              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1112  #pragma omp for nowait  #pragma omp for nowait
1113                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 725  void Brick::setToGradient(escript::Data& Line 1125  void Brick::setToGradient(escript::Data&
1125                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / h0;
1126                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;
1127                              o[INDEX3(i,2,0,numComp,3)] = (f_011[i]+f_111[i]-f_010[i]-f_110[i])*C4 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_011[i]+f_111[i]-f_010[i]-f_110[i])*C4 / h2;
1128                          } /* end of component loop i */                          } // end of component loop i
1129                      } /* end of k0 loop */                      } // end of k0 loop
1130                  } /* end of k2 loop */                  } // end of k2 loop
1131              } /* end of face 3 */              } // end of face 3
1132              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1133  #pragma omp for nowait  #pragma omp for nowait
1134                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 746  void Brick::setToGradient(escript::Data& Line 1146  void Brick::setToGradient(escript::Data&
1146                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / h0;
1147                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / h1;
1148                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C4 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C4 / h2;
1149                          } /* end of component loop i */                          } // end of component loop i
1150                      } /* end of k0 loop */                      } // end of k0 loop
1151                  } /* end of k1 loop */                  } // end of k1 loop
1152              } /* end of face 4 */              } // end of face 4
1153              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1154  #pragma omp for nowait  #pragma omp for nowait
1155                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 767  void Brick::setToGradient(escript::Data& Line 1167  void Brick::setToGradient(escript::Data&
1167                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / h0;
1168                              o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / h1;                              o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / h1;
1169                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / h2;                              o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / h2;
1170                          } /* end of component loop i */                          } // end of component loop i
1171                      } /* end of k0 loop */                      } // end of k0 loop
1172                  } /* end of k1 loop */                  } // end of k1 loop
1173              } /* end of face 5 */              } // end of face 5
1174          } // end of parallel section          } // end of parallel section
         /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */  
     } else {  
         stringstream msg;  
         msg << "setToGradient() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
1175      }      }
1176  }  }
1177    
1178  void Brick::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  //protected
1179    void Brick::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1180  {  {
1181      escript::Data& in = *const_cast<escript::Data*>(&arg);      const dim_t numComp = arg.getDataPointSize();
     const dim_t numComp = in.getDataPointSize();  
1182      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1183      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
1184      const double h2 = m_l2/m_gNE2;      const double h2 = m_l2/m_gNE2;
1185        const index_t left = (m_offset0==0 ? 0 : 1);
1186        const index_t bottom = (m_offset1==0 ? 0 : 1);
1187        const index_t front = (m_offset2==0 ? 0 : 1);
1188      if (arg.getFunctionSpace().getTypeCode() == Elements) {      if (arg.getFunctionSpace().getTypeCode() == Elements) {
1189          const double w_0 = h0*h1*h2/8.;          const double w_0 = h0*h1*h2/8.;
1190  #pragma omp parallel  #pragma omp parallel
1191          {          {
1192              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1193  #pragma omp for nowait  #pragma omp for nowait
1194              for (index_t k2 = 0; k2 < m_NE2; ++k2) {              for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1195                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1196                      for (index_t k0 = 0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1197                          const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));                          const double* f = arg.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
1198                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1199                              const register double f_0 = f[INDEX2(i,0,numComp)];                              const register double f_0 = f[INDEX2(i,0,numComp)];
1200                              const register double f_1 = f[INDEX2(i,1,numComp)];                              const register double f_1 = f[INDEX2(i,1,numComp)];
# Line 808  void Brick::setToIntegrals(vector<double Line 1205  void Brick::setToIntegrals(vector<double
1205                              const register double f_6 = f[INDEX2(i,6,numComp)];                              const register double f_6 = f[INDEX2(i,6,numComp)];
1206                              const register double f_7 = f[INDEX2(i,7,numComp)];                              const register double f_7 = f[INDEX2(i,7,numComp)];
1207                              int_local[i]+=(f_0+f_1+f_2+f_3+f_4+f_5+f_6+f_7)*w_0;                              int_local[i]+=(f_0+f_1+f_2+f_3+f_4+f_5+f_6+f_7)*w_0;
1208                          }  /* end of component loop i */                          }  // end of component loop i
1209                      } /* end of k0 loop */                      } // end of k0 loop
1210                  } /* end of k1 loop */                  } // end of k1 loop
1211              } /* end of k2 loop */              } // end of k2 loop
1212    
1213  #pragma omp critical  #pragma omp critical
1214              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
# Line 823  void Brick::setToIntegrals(vector<double Line 1220  void Brick::setToIntegrals(vector<double
1220          {          {
1221              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1222  #pragma omp for nowait  #pragma omp for nowait
1223              for (index_t k2 = 0; k2 < m_NE2; ++k2) {              for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1224                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1225                      for (index_t k0 = 0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1226                          const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));                          const double* f = arg.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
1227                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1228                              int_local[i]+=f[i]*w_0;                              int_local[i]+=f[i]*w_0;
1229                          }  /* end of component loop i */                          }  // end of component loop i
1230                      } /* end of k0 loop */                      } // end of k0 loop
1231                  } /* end of k1 loop */                  } // end of k1 loop
1232              } /* end of k2 loop */              } // end of k2 loop
1233    
1234  #pragma omp critical  #pragma omp critical
1235              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
# Line 847  void Brick::setToIntegrals(vector<double Line 1244  void Brick::setToIntegrals(vector<double
1244              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1245              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1246  #pragma omp for nowait  #pragma omp for nowait
1247                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1248                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1249                          const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          const double* f = arg.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1250                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1251                              const register double f_0 = f[INDEX2(i,0,numComp)];                              const register double f_0 = f[INDEX2(i,0,numComp)];
1252                              const register double f_1 = f[INDEX2(i,1,numComp)];                              const register double f_1 = f[INDEX2(i,1,numComp)];
1253                              const register double f_2 = f[INDEX2(i,2,numComp)];                              const register double f_2 = f[INDEX2(i,2,numComp)];
1254                              const register double f_3 = f[INDEX2(i,3,numComp)];                              const register double f_3 = f[INDEX2(i,3,numComp)];
1255                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
1256                          }  /* end of component loop i */                          }  // end of component loop i
1257                      } /* end of k1 loop */                      } // end of k1 loop
1258                  } /* end of k2 loop */                  } // end of k2 loop
1259              }              }
1260    
1261              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1262  #pragma omp for nowait  #pragma omp for nowait
1263                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1264                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1265                          const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          const double* f = arg.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1266                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1267                              const register double f_0 = f[INDEX2(i,0,numComp)];                              const register double f_0 = f[INDEX2(i,0,numComp)];
1268                              const register double f_1 = f[INDEX2(i,1,numComp)];                              const register double f_1 = f[INDEX2(i,1,numComp)];
1269                              const register double f_2 = f[INDEX2(i,2,numComp)];                              const register double f_2 = f[INDEX2(i,2,numComp)];
1270                              const register double f_3 = f[INDEX2(i,3,numComp)];                              const register double f_3 = f[INDEX2(i,3,numComp)];
1271                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
1272                          }  /* end of component loop i */                          }  // end of component loop i
1273                      } /* end of k1 loop */                      } // end of k1 loop
1274                  } /* end of k2 loop */                  } // end of k2 loop
1275              }              }
1276    
1277              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1278  #pragma omp for nowait  #pragma omp for nowait
1279                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1280                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1281                          const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          const double* f = arg.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1282                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1283                              const register double f_0 = f[INDEX2(i,0,numComp)];                              const register double f_0 = f[INDEX2(i,0,numComp)];
1284                              const register double f_1 = f[INDEX2(i,1,numComp)];                              const register double f_1 = f[INDEX2(i,1,numComp)];
1285                              const register double f_2 = f[INDEX2(i,2,numComp)];                              const register double f_2 = f[INDEX2(i,2,numComp)];
1286                              const register double f_3 = f[INDEX2(i,3,numComp)];                              const register double f_3 = f[INDEX2(i,3,numComp)];
1287                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
1288                          }  /* end of component loop i */                          }  // end of component loop i
1289                      } /* end of k1 loop */                      } // end of k1 loop
1290                  } /* end of k2 loop */                  } // end of k2 loop
1291              }              }
1292    
1293              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1294  #pragma omp for nowait  #pragma omp for nowait
1295                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1296                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1297                          const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          const double* f = arg.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1298                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1299                              const register double f_0 = f[INDEX2(i,0,numComp)];                              const register double f_0 = f[INDEX2(i,0,numComp)];
1300                              const register double f_1 = f[INDEX2(i,1,numComp)];                              const register double f_1 = f[INDEX2(i,1,numComp)];
1301                              const register double f_2 = f[INDEX2(i,2,numComp)];                              const register double f_2 = f[INDEX2(i,2,numComp)];
1302                              const register double f_3 = f[INDEX2(i,3,numComp)];                              const register double f_3 = f[INDEX2(i,3,numComp)];
1303                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
1304                          }  /* end of component loop i */                          }  // end of component loop i
1305                      } /* end of k1 loop */                      } // end of k1 loop
1306                  } /* end of k2 loop */                  } // end of k2 loop
1307              }              }
1308    
1309              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1310  #pragma omp for nowait  #pragma omp for nowait
1311                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1312                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1313                          const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          const double* f = arg.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1314                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1315                              const register double f_0 = f[INDEX2(i,0,numComp)];                              const register double f_0 = f[INDEX2(i,0,numComp)];
1316                              const register double f_1 = f[INDEX2(i,1,numComp)];                              const register double f_1 = f[INDEX2(i,1,numComp)];
1317                              const register double f_2 = f[INDEX2(i,2,numComp)];                              const register double f_2 = f[INDEX2(i,2,numComp)];
1318                              const register double f_3 = f[INDEX2(i,3,numComp)];                              const register double f_3 = f[INDEX2(i,3,numComp)];
1319                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
1320                          }  /* end of component loop i */                          }  // end of component loop i
1321                      } /* end of k1 loop */                      } // end of k1 loop
1322                  } /* end of k2 loop */                  } // end of k2 loop
1323              }              }
1324    
1325              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1326  #pragma omp for nowait  #pragma omp for nowait
1327                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1328                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1329                          const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          const double* f = arg.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1330                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1331                              const register double f_0 = f[INDEX2(i,0,numComp)];                              const register double f_0 = f[INDEX2(i,0,numComp)];
1332                              const register double f_1 = f[INDEX2(i,1,numComp)];                              const register double f_1 = f[INDEX2(i,1,numComp)];
1333                              const register double f_2 = f[INDEX2(i,2,numComp)];                              const register double f_2 = f[INDEX2(i,2,numComp)];
1334                              const register double f_3 = f[INDEX2(i,3,numComp)];                              const register double f_3 = f[INDEX2(i,3,numComp)];
1335                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;                              int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
1336                          }  /* end of component loop i */                          }  // end of component loop i
1337                      } /* end of k1 loop */                      } // end of k1 loop
1338                  } /* end of k2 loop */                  } // end of k2 loop
1339              }              }
1340    
1341  #pragma omp critical  #pragma omp critical
# Line 955  void Brick::setToIntegrals(vector<double Line 1352  void Brick::setToIntegrals(vector<double
1352              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1353              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1354  #pragma omp for nowait  #pragma omp for nowait
1355                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1356                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1357                          const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          const double* f = arg.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1358                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1359                              int_local[i]+=f[i]*w_0;                              int_local[i]+=f[i]*w_0;
1360                          }  /* end of component loop i */                          }  // end of component loop i
1361                      } /* end of k1 loop */                      } // end of k1 loop
1362                  } /* end of k2 loop */                  } // end of k2 loop
1363              }              }
1364    
1365              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1366  #pragma omp for nowait  #pragma omp for nowait
1367                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1368                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1369                          const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          const double* f = arg.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1370                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1371                              int_local[i]+=f[i]*w_0;                              int_local[i]+=f[i]*w_0;
1372                          }  /* end of component loop i */                          }  // end of component loop i
1373                      } /* end of k1 loop */                      } // end of k1 loop
1374                  } /* end of k2 loop */                  } // end of k2 loop
1375              }              }
1376    
1377              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1378  #pragma omp for nowait  #pragma omp for nowait
1379                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1380                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1381                          const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          const double* f = arg.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1382                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1383                              int_local[i]+=f[i]*w_1;                              int_local[i]+=f[i]*w_1;
1384                          }  /* end of component loop i */                          }  // end of component loop i
1385                      } /* end of k1 loop */                      } // end of k1 loop
1386                  } /* end of k2 loop */                  } // end of k2 loop
1387              }              }
1388    
1389              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1390  #pragma omp for nowait  #pragma omp for nowait
1391                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2 = front; k2 < front+m_ownNE2; ++k2) {
1392                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1393                          const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          const double* f = arg.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1394                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1395                              int_local[i]+=f[i]*w_1;                              int_local[i]+=f[i]*w_1;
1396                          }  /* end of component loop i */                          }  // end of component loop i
1397                      } /* end of k1 loop */                      } // end of k1 loop
1398                  } /* end of k2 loop */                  } // end of k2 loop
1399              }              }
1400    
1401              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1402  #pragma omp for nowait  #pragma omp for nowait
1403                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1404                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1405                          const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          const double* f = arg.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1406                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1407                              int_local[i]+=f[i]*w_2;                              int_local[i]+=f[i]*w_2;
1408                          }  /* end of component loop i */                          }  // end of component loop i
1409                      } /* end of k1 loop */                      } // end of k1 loop
1410                  } /* end of k2 loop */                  } // end of k2 loop
1411              }              }
1412    
1413              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1414  #pragma omp for nowait  #pragma omp for nowait
1415                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
1416                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
1417                          const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          const double* f = arg.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1418                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1419                              int_local[i]+=f[i]*w_2;                              int_local[i]+=f[i]*w_2;
1420                          }  /* end of component loop i */                          }  // end of component loop i
1421                      } /* end of k1 loop */                      } // end of k1 loop
1422                  } /* end of k2 loop */                  } // end of k2 loop
1423              }              }
1424    
1425  #pragma omp critical  #pragma omp critical
# Line 1030  void Brick::setToIntegrals(vector<double Line 1427  void Brick::setToIntegrals(vector<double
1427                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1428          } // end of parallel section          } // end of parallel section
1429    
     } else {  
         stringstream msg;  
         msg << "setToIntegrals() not implemented for "  
             << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
 }  
   
 void Brick::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 k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                         double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));  
                         // set vector at four quadrature points  
                         *o++ = -1.; *o++ = 0.; *o++ = 0.;  
                         *o++ = -1.; *o++ = 0.; *o++ = 0.;  
                         *o++ = -1.; *o++ = 0.; *o++ = 0.;  
                         *o++ = -1.; *o++ = 0.; *o = 0.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                         double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));  
                         // set vector at four quadrature points  
                         *o++ = 1.; *o++ = 0.; *o++ = 0.;  
                         *o++ = 1.; *o++ = 0.; *o++ = 0.;  
                         *o++ = 1.; *o++ = 0.; *o++ = 0.;  
                         *o++ = 1.; *o++ = 0.; *o = 0.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));  
                         // set vector at four quadrature points  
                         *o++ = 0.; *o++ = -1.; *o++ = 0.;  
                         *o++ = 0.; *o++ = -1.; *o++ = 0.;  
                         *o++ = 0.; *o++ = -1.; *o++ = 0.;  
                         *o++ = 0.; *o++ = -1.; *o = 0.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));  
                         // set vector at four quadrature points  
                         *o++ = 0.; *o++ = 1.; *o++ = 0.;  
                         *o++ = 0.; *o++ = 1.; *o++ = 0.;  
                         *o++ = 0.; *o++ = 1.; *o++ = 0.;  
                         *o++ = 0.; *o++ = 1.; *o = 0.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[4] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));  
                         // set vector at four quadrature points  
                         *o++ = 0.; *o++ = 0.; *o++ = -1.;  
                         *o++ = 0.; *o++ = 0.; *o++ = -1.;  
                         *o++ = 0.; *o++ = 0.; *o++ = -1.;  
                         *o++ = 0.; *o++ = 0.; *o = -1.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[5] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));  
                         // set vector at four quadrature points  
                         *o++ = 0.; *o++ = 0.; *o++ = 1.;  
                         *o++ = 0.; *o++ = 0.; *o++ = 1.;  
                         *o++ = 0.; *o++ = 0.; *o++ = 1.;  
                         *o++ = 0.; *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 k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                         double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));  
                         *o++ = -1.;  
                         *o++ = 0.;  
                         *o = 0.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[1] > -1) {  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                         double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));  
                         *o++ = 1.;  
                         *o++ = 0.;  
                         *o = 0.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[2] > -1) {  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));  
                         *o++ = 0.;  
                         *o++ = -1.;  
                         *o = 0.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[3] > -1) {  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));  
                         *o++ = 0.;  
                         *o++ = 1.;  
                         *o = 0.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[4] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));  
                         *o++ = 0.;  
                         *o++ = 0.;  
                         *o = -1.;  
                     }  
                 }  
             }  
   
             if (m_faceOffset[5] > -1) {  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));  
                         *o++ = 0.;  
                         *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 Brick::setToSize(escript::Data& out) const  
 {  
     if (out.getFunctionSpace().getTypeCode() == Elements  
             || out.getFunctionSpace().getTypeCode() == ReducedElements) {  
         out.requireWrite();  
         const dim_t numQuad=out.getNumDataPointsPerSample();  
         const double xSize=getFirstCoordAndSpacing(0).second;  
         const double ySize=getFirstCoordAndSpacing(1).second;  
         const double zSize=getFirstCoordAndSpacing(2).second;  
         const double size=min(min(xSize,ySize),zSize);  
 #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 xSize=getFirstCoordAndSpacing(0).second;  
         const double ySize=getFirstCoordAndSpacing(1).second;  
         const double zSize=getFirstCoordAndSpacing(2).second;  
 #pragma omp parallel  
         {  
             if (m_faceOffset[0] > -1) {  
                 const double size=min(ySize,zSize);  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                         double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));  
                         fill(o, o+numQuad, size);  
                     }  
                 }  
             }  
   
             if (m_faceOffset[1] > -1) {  
                 const double size=min(ySize,zSize);  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                         double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));  
                         fill(o, o+numQuad, size);  
                     }  
                 }  
             }  
   
             if (m_faceOffset[2] > -1) {  
                 const double size=min(xSize,zSize);  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));  
                         fill(o, o+numQuad, size);  
                     }  
                 }  
             }  
   
             if (m_faceOffset[3] > -1) {  
                 const double size=min(xSize,zSize);  
 #pragma omp for nowait  
                 for (index_t k2 = 0; k2 < m_NE2; ++k2) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));  
                         fill(o, o+numQuad, size);  
                     }  
                 }  
             }  
   
             if (m_faceOffset[4] > -1) {  
                 const double size=min(xSize,ySize);  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));  
                         fill(o, o+numQuad, size);  
                     }  
                 }  
             }  
   
             if (m_faceOffset[5] > -1) {  
                 const double size=min(xSize,ySize);  
 #pragma omp for nowait  
                 for (index_t k1 = 0; k1 < m_NE1; ++k1) {  
                     for (index_t k0 = 0; k0 < m_NE0; ++k0) {  
                         double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));  
                         fill(o, o+numQuad, size);  
                     }  
                 }  
             }  
         } // end of parallel section  
   
     } else {  
         stringstream msg;  
         msg << "setToSize() not implemented for "  
             << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());  
         throw RipleyException(msg.str());  
     }  
 }  
   
 Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,  
                                             bool reducedColOrder) const  
 {  
     if (reducedRowOrder || reducedColOrder)  
         throw RipleyException("getPattern() not implemented for reduced order");  
   
     return m_pattern;  
 }  
   
 void Brick::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);  
         pair<double,double> zdz = getFirstCoordAndSpacing(2);  
         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*m_N1)/m_N0)*ydy.second  
                 << "  " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;  
         }  
     }  
 }  
   
 IndexVector Brick::getNumNodesPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_N0);  
     ret.push_back(m_N1);  
     ret.push_back(m_N2);  
     return ret;  
 }  
   
 IndexVector Brick::getNumElementsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NE0);  
     ret.push_back(m_NE1);  
     ret.push_back(m_NE2);  
     return ret;  
 }  
   
 IndexVector Brick::getNumFacesPerBoundary() const  
 {  
     IndexVector ret(6, 0);  
     //left  
     if (m_offset0==0)  
         ret[0]=m_NE1*m_NE2;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         ret[1]=m_NE1*m_NE2;  
     //bottom  
     if (m_offset1==0)  
         ret[2]=m_NE0*m_NE2;  
     //top  
     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)  
         ret[3]=m_NE0*m_NE2;  
     //front  
     if (m_offset2==0)  
         ret[4]=m_NE0*m_NE1;  
     //back  
     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)  
         ret[5]=m_NE0*m_NE1;  
     return ret;  
 }  
   
 pair<double,double> Brick::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);  
     else if (dim==2)  
         return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);  
   
     throw RipleyException("getFirstCoordAndSpacing(): invalid argument");  
 }  
   
 //protected  
 dim_t Brick::getNumDOF() const  
 {  
     return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY*(m_gNE2+1)/m_NZ;  
 }  
   
 //protected  
 dim_t Brick::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 Brick::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);  
     pair<double,double> zdz = getFirstCoordAndSpacing(2);  
     arg.requireWrite();  
 #pragma omp parallel for  
     for (dim_t i2 = 0; i2 < m_N2; i2++) {  
         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+m_N0*m_N1*i2);  
                 point[0] = xdx.first+i0*xdx.second;  
                 point[1] = ydy.first+i1*ydy.second;  
                 point[2] = zdz.first+i2*zdz.second;  
             }  
         }  
1430      }      }
1431  }  }
1432    
# Line 1918  void Brick::interpolateNodesOnElements(e Line 1911  void Brick::interpolateNodesOnElements(e
1911      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1912      if (reduced) {      if (reduced) {
1913          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */  
1914          const double c0 = .125;          const double c0 = .125;
1915  #pragma omp parallel for  #pragma omp parallel for
1916          for (index_t k2=0; k2 < m_NE2; ++k2) {          for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 1935  void Brick::interpolateNodesOnElements(e Line 1927  void Brick::interpolateNodesOnElements(e
1927                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1928                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1929                          o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i] + f_100[i] + f_101[i] + f_110[i] + f_111[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i] + f_100[i] + f_101[i] + f_110[i] + f_111[i]);
1930                      } /* end of component loop i */                      } // end of component loop i
1931                  } /* end of k0 loop */                  } // end of k0 loop
1932              } /* end of k1 loop */              } // end of k1 loop
1933          } /* end of k2 loop */          } // end of k2 loop
         /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */  
1934      } else {      } else {
1935          out.requireWrite();          out.requireWrite();
         /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */  
1936          const double c0 = .0094373878376559314545;          const double c0 = .0094373878376559314545;
1937          const double c1 = .035220810900864519624;          const double c1 = .035220810900864519624;
1938          const double c2 = .13144585576580214704;          const double c2 = .13144585576580214704;
# Line 1969  void Brick::interpolateNodesOnElements(e Line 1959  void Brick::interpolateNodesOnElements(e
1959                          o[INDEX2(i,numComp,5)] = f_010[i]*c0 + f_101[i]*c3 + c2*(f_001[i] + f_100[i] + f_111[i]) + c1*(f_000[i] + f_011[i] + f_110[i]);                          o[INDEX2(i,numComp,5)] = f_010[i]*c0 + f_101[i]*c3 + c2*(f_001[i] + f_100[i] + f_111[i]) + c1*(f_000[i] + f_011[i] + f_110[i]);
1960                          o[INDEX2(i,numComp,6)] = f_011[i]*c3 + f_100[i]*c0 + c2*(f_001[i] + f_010[i] + f_111[i]) + c1*(f_000[i] + f_101[i] + f_110[i]);                          o[INDEX2(i,numComp,6)] = f_011[i]*c3 + f_100[i]*c0 + c2*(f_001[i] + f_010[i] + f_111[i]) + c1*(f_000[i] + f_101[i] + f_110[i]);
1961                          o[INDEX2(i,numComp,7)] = f_000[i]*c0 + f_111[i]*c3 + c2*(f_011[i] + f_101[i] + f_110[i]) + c1*(f_001[i] + f_010[i] + f_100[i]);                          o[INDEX2(i,numComp,7)] = f_000[i]*c0 + f_111[i]*c3 + c2*(f_011[i] + f_101[i] + f_110[i]) + c1*(f_001[i] + f_010[i] + f_100[i]);
1962                      } /* end of component loop i */                      } // end of component loop i
1963                  } /* end of k0 loop */                  } // end of k0 loop
1964              } /* end of k1 loop */              } // end of k1 loop
1965          } /* end of k2 loop */          } // end of k2 loop
         /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */  
1966      }      }
1967  }  }
1968    
# Line 1987  void Brick::interpolateNodesOnFaces(escr Line 1976  void Brick::interpolateNodesOnFaces(escr
1976          const double c0 = .25;          const double c0 = .25;
1977  #pragma omp parallel  #pragma omp parallel
1978          {          {
             /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */  
1979              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1980  #pragma omp for nowait  #pragma omp for nowait
1981                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 1999  void Brick::interpolateNodesOnFaces(escr Line 1987  void Brick::interpolateNodesOnFaces(escr
1987                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1988                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1989                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]);
1990                          } /* end of component loop i */                          } // end of component loop i
1991                      } /* end of k1 loop */                      } // end of k1 loop
1992                  } /* end of k2 loop */                  } // end of k2 loop
1993              } /* end of face 0 */              } // end of face 0
1994              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1995  #pragma omp for nowait  #pragma omp for nowait
1996                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 2014  void Brick::interpolateNodesOnFaces(escr Line 2002  void Brick::interpolateNodesOnFaces(escr
2002                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
2003                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2004                              o[INDEX2(i,numComp,0)] = c0*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
2005                          } /* end of component loop i */                          } // end of component loop i
2006                      } /* end of k1 loop */                      } // end of k1 loop
2007                  } /* end of k2 loop */                  } // end of k2 loop
2008              } /* end of face 1 */              } // end of face 1
2009              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
2010  #pragma omp for nowait  #pragma omp for nowait
2011                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 2029  void Brick::interpolateNodesOnFaces(escr Line 2017  void Brick::interpolateNodesOnFaces(escr
2017                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
2018                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2019                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);
2020                          } /* end of component loop i */                          } // end of component loop i
2021                      } /* end of k0 loop */                      } // end of k0 loop
2022                  } /* end of k2 loop */                  } // end of k2 loop
2023              } /* end of face 2 */              } // end of face 2
2024              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
2025  #pragma omp for nowait  #pragma omp for nowait
2026                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 2044  void Brick::interpolateNodesOnFaces(escr Line 2032  void Brick::interpolateNodesOnFaces(escr
2032                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
2033                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2034                              o[INDEX2(i,numComp,0)] = c0*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);
2035                          } /* end of component loop i */                          } // end of component loop i
2036                      } /* end of k0 loop */                      } // end of k0 loop
2037                  } /* end of k2 loop */                  } // end of k2 loop
2038              } /* end of face 3 */              } // end of face 3
2039              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
2040  #pragma omp for nowait  #pragma omp for nowait
2041                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 2059  void Brick::interpolateNodesOnFaces(escr Line 2047  void Brick::interpolateNodesOnFaces(escr
2047                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
2048                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2049                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);
2050                          } /* end of component loop i */                          } // end of component loop i
2051                      } /* end of k0 loop */                      } // end of k0 loop
2052                  } /* end of k1 loop */                  } // end of k1 loop
2053              } /* end of face 4 */              } // end of face 4
2054              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
2055  #pragma omp for nowait  #pragma omp for nowait
2056                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 2074  void Brick::interpolateNodesOnFaces(escr Line 2062  void Brick::interpolateNodesOnFaces(escr
2062                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
2063                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2064                              o[INDEX2(i,numComp,0)] = c0*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);
2065                          } /* end of component loop i */                          } // end of component loop i
2066                      } /* end of k0 loop */                      } // end of k0 loop
2067                  } /* end of k1 loop */                  } // end of k1 loop
2068              } /* end of face 5 */              } // end of face 5
             /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */  
2069          } // end of parallel section          } // end of parallel section
2070      } else {      } else {
2071          out.requireWrite();          out.requireWrite();
# Line 2087  void Brick::interpolateNodesOnFaces(escr Line 2074  void Brick::interpolateNodesOnFaces(escr
2074          const double c2 = 0.62200846792814621559;          const double c2 = 0.62200846792814621559;
2075  #pragma omp parallel  #pragma omp parallel
2076          {          {
             /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */  
2077              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
2078  #pragma omp for nowait  #pragma omp for nowait
2079                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 2102  void Brick::interpolateNodesOnFaces(escr Line 2088  void Brick::interpolateNodesOnFaces(escr
2088                              o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_010[i]*c2 + c1*(f_000[i] + f_011[i]);                              o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_010[i]*c2 + c1*(f_000[i] + f_011[i]);
2089                              o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_010[i]*c0 + c1*(f_000[i] + f_011[i]);                              o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_010[i]*c0 + c1*(f_000[i] + f_011[i]);
2090                              o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_011[i]*c2 + c1*(f_001[i] + f_010[i]);                              o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_011[i]*c2 + c1*(f_001[i] + f_010[i]);
2091                          } /* end of component loop i */                          } // end of component loop i
2092                      } /* end of k1 loop */                      } // end of k1 loop
2093                  } /* end of k2 loop */                  } // end of k2 loop
2094              } /* end of face 0 */              } // end of face 0
2095              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
2096  #pragma omp for nowait  #pragma omp for nowait
2097                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 2120  void Brick::interpolateNodesOnFaces(escr Line 2106  void Brick::interpolateNodesOnFaces(escr
2106                              o[INDEX2(i,numComp,1)] = f_101[i]*c0 + f_110[i]*c2 + c1*(f_100[i] + f_111[i]);                              o[INDEX2(i,numComp,1)] = f_101[i]*c0 + f_110[i]*c2 + c1*(f_100[i] + f_111[i]);
2107                              o[INDEX2(i,numComp,2)] = f_101[i]*c2 + f_110[i]*c0 + c1*(f_100[i] + f_111[i]);                              o[INDEX2(i,numComp,2)] = f_101[i]*c2 + f_110[i]*c0 + c1*(f_100[i] + f_111[i]);
2108                              o[INDEX2(i,numComp,3)] = f_100[i]*c0 + f_111[i]*c2 + c1*(f_101[i] + f_110[i]);                              o[INDEX2(i,numComp,3)] = f_100[i]*c0 + f_111[i]*c2 + c1*(f_101[i] + f_110[i]);
2109                          } /* end of component loop i */                          } // end of component loop i
2110                      } /* end of k1 loop */                      } // end of k1 loop
2111                  } /* end of k2 loop */                  } // end of k2 loop
2112              } /* end of face 1 */              } // end of face 1
2113              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
2114  #pragma omp for nowait  #pragma omp for nowait
2115                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 2138  void Brick::interpolateNodesOnFaces(escr Line 2124  void Brick::interpolateNodesOnFaces(escr
2124                              o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_101[i]);                              o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_101[i]);
2125                              o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_101[i]);                              o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_101[i]);
2126                              o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_100[i]);                              o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_100[i]);
2127                          } /* end of component loop i */                          } // end of component loop i
2128                      } /* end of k0 loop */                      } // end of k0 loop
2129                  } /* end of k2 loop */                  } // end of k2 loop
2130              } /* end of face 2 */              } // end of face 2
2131              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
2132  #pragma omp for nowait  #pragma omp for nowait
2133                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
# Line 2156  void Brick::interpolateNodesOnFaces(escr Line 2142  void Brick::interpolateNodesOnFaces(escr
2142                              o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_111[i]);                              o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_111[i]);
2143                              o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_111[i]);                              o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_111[i]);
2144                              o[INDEX2(i,numComp,3)] = f_010[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_110[i]);                              o[INDEX2(i,numComp,3)] = f_010[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_110[i]);
2145                          } /* end of component loop i */                          } // end of component loop i
2146                      } /* end of k0 loop */                      } // end of k0 loop
2147                  } /* end of k2 loop */                  } // end of k2 loop
2148              } /* end of face 3 */              } // end of face 3
2149              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
2150  #pragma omp for nowait  #pragma omp for nowait
2151                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 2174  void Brick::interpolateNodesOnFaces(escr Line 2160  void Brick::interpolateNodesOnFaces(escr
2160                              o[INDEX2(i,numComp,1)] = f_010[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_110[i]);                              o[INDEX2(i,numComp,1)] = f_010[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_110[i]);
2161                              o[INDEX2(i,numComp,2)] = f_010[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_110[i]);                              o[INDEX2(i,numComp,2)] = f_010[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_110[i]);
2162                              o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_100[i]);                              o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_100[i]);
2163                          } /* end of component loop i */                          } // end of component loop i
2164                      } /* end of k0 loop */                      } // end of k0 loop
2165                  } /* end of k1 loop */                  } // end of k1 loop
2166              } /* end of face 4 */              } // end of face 4
2167              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
2168  #pragma omp for nowait  #pragma omp for nowait
2169                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 2192  void Brick::interpolateNodesOnFaces(escr Line 2178  void Brick::interpolateNodesOnFaces(escr
2178                              o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_111[i]);                              o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_111[i]);
2179                              o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_111[i]);                              o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_111[i]);
2180                              o[INDEX2(i,numComp,3)] = f_001[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_101[i]);                              o[INDEX2(i,numComp,3)] = f_001[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_101[i]);
2181                          } /* end of component loop i */                          } // end of component loop i
2182                      } /* end of k0 loop */                      } // end of k0 loop
2183                  } /* end of k1 loop */                  } // end of k1 loop
2184              } /* end of face 5 */              } // end of face 5
             /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */  
2185          } // end of parallel section          } // end of parallel section
2186      }      }
2187  }  }
# Line 2211  void Brick::assemblePDESingle(Paso_Syste Line 2196  void Brick::assemblePDESingle(Paso_Syste
2196      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
2197      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
2198      const double h2 = m_l2/m_gNE2;      const double h2 = m_l2/m_gNE2;
     /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */  
2199      const double w0 = 0.0009303791403858427308*h1*h2/h0;      const double w0 = 0.0009303791403858427308*h1*h2/h0;
2200      const double w1 = 0.0009303791403858427308*h2;      const double w1 = 0.0009303791403858427308*h2;
2201        const double w2 = -0.00024929433932114870101*h1;
2202        const double w3 = 0.0009303791403858427308*h0*h2/h1;
2203        const double w4 = -0.00024929433932114870101*h0;
2204        const double w5 = 0.0009303791403858427308*h1;
2205        const double w6 = 0.0009303791403858427308*h0;
2206        const double w7 = -0.00024929433932114870101*h0*h1/h2;
2207        const double w8 = 0.0034722222222222222222*h2;
2208        const double w9 = -0.0009303791403858427308*h1;
2209      const double w10 = 0.012958509748503046158*h0*h2/h1;      const double w10 = 0.012958509748503046158*h0*h2/h1;
     const double w100 = 0.038141762351741127649*h0*h1;  
     const double w101 = 0.000052682092703316795705*h0*h1;  
     const double w102 = 0.0007337668937680108255*h0*h1;  
     const double w103 = 0.010220054420048834761*h0*h1;  
     const double w104 = -0.0001966122466178319053*h1*h2;  
     const double w105 = -0.0001966122466178319053*h0*h2;  
     const double w106 = -0.0007337668937680108255*h1*h2;  
     const double w107 = -0.0007337668937680108255*h0*h2;  
     const double w108 = -0.0027384553284542113967*h1*h2;  
     const double w109 = -0.0027384553284542113967*h0*h2;  
2210      const double w11 = -0.0034722222222222222222*h0;      const double w11 = -0.0034722222222222222222*h0;
     const double w110 = -0.010220054420048834761*h1*h2;  
     const double w111 = -0.010220054420048834761*h0*h2;  
     const double w112 = -0.0007337668937680108255*h0*h1;  
     const double w113 = -0.010220054420048834761*h0*h1;  
     const double w114 = -0.038141762351741127649*h0*h2;  
     const double w115 = -0.000052682092703316795705*h0*h2;  
     const double w116 = -0.0001966122466178319053*h0*h1;  
     const double w117 = -0.0027384553284542113967*h0*h1;  
     const double w118 = 0.000052682092703316795705*h0*h2;  
     const double w119 = 0.038141762351741127649*h0*h2;  
2211      const double w12 = 0.0034722222222222222222*h1;      const double w12 = 0.0034722222222222222222*h1;
     const double w120 = 0.000052682092703316795705*h1*h2;  
     const double w121 = 0.038141762351741127649*h1*h2;  
     const double w122 = -0.000052682092703316795705*h0*h1;  
     const double w123 = -0.038141762351741127649*h0*h1;  
     const double w124 = -0.000052682092703316795705*h1*h2;  
     const double w125 = -0.038141762351741127649*h1*h2;  
     const double w126 = 0.027777777777777777778*h1*h2;  
     const double w127 = 0.027777777777777777778*h0*h2;  
     const double w128 = 0.055555555555555555556*h0*h1;  
     const double w129 = -0.027777777777777777778*h1*h2;  
2212      const double w13 = 0.012958509748503046158*h0;      const double w13 = 0.012958509748503046158*h0;
     const double w130 = -0.027777777777777777778*h0*h2;  
     const double w131 = 0.013888888888888888889*h0*h1;  
     const double w132 = -0.055555555555555555556*h0*h2;  
     const double w133 = -0.027777777777777777778*h0*h1;  
     const double w134 = 0.055555555555555555556*h0*h2;  
     const double w135 = 0.027777777777777777778*h0*h1;  
     const double w136 = -0.013888888888888888889*h0*h1;  
     const double w137 = 0.055555555555555555556*h1*h2;  
     const double w138 = -0.013888888888888888889*h1*h2;  
     const double w139 = -0.013888888888888888889*h0*h2;  
2213      const double w14 = -0.0034722222222222222222*h0*h1/h2;      const double w14 = -0.0034722222222222222222*h0*h1/h2;
     const double w140 = -0.055555555555555555556*h0*h1;  
     const double w141 = 0.013888888888888888889*h1*h2;  
     const double w142 = 0.013888888888888888889*h0*h2;  
     const double w143 = -0.055555555555555555556*h1*h2;  
     const double w144 = 0.000041549056553524783501*h0*h1*h2;  
     const double w145 = 0.0005787037037037037037*h0*h1*h2;  
     const double w146 = 0.0080603027952983270684*h0*h1*h2;  
     const double w147 = 0.0001550631900643071218*h0*h1*h2;  
     const double w148 = 0.002159751624750507693*h0*h1*h2;  
     const double w149 = 0.03008145955644280058*h0*h1*h2;  
2214      const double w15 = 0.012958509748503046158*h1*h2/h0;      const double w15 = 0.012958509748503046158*h1*h2/h0;
     const double w150 = 0.000011133036149792012204*h0*h1*h2;  
     const double w151 = 0.018518518518518518519*h0*h1*h2;  
     const double w152 = 0.0092592592592592592592*h0*h1*h2;  
     const double w153 = 0.0046296296296296296296*h0*h1*h2;  
     const double w154 = 0.037037037037037037037*h0*h1*h2;  
     const double w155 = -0.077751058491018276949*h1*h2;  
     const double w156 = -0.077751058491018276949*h0*h2;  
     const double w157 = -0.077751058491018276949*h0*h1;  
     const double w158 = -0.020833333333333333333*h0*h2;  
     const double w159 = -0.020833333333333333333*h0*h1;  
2215      const double w16 = -0.0034722222222222222222*h1;      const double w16 = -0.0034722222222222222222*h1;
     const double w160 = -0.020833333333333333333*h1*h2;  
     const double w161 = -0.0055822748423150563848*h0*h1;  
     const double w162 = -0.0055822748423150563848*h0*h2;  
     const double w163 = -0.0055822748423150563848*h1*h2;  
     const double w164 = 0.077751058491018276949*h1*h2;  
     const double w165 = 0.020833333333333333333*h1*h2;  
     const double w166 = 0.0055822748423150563848*h1*h2;  
     const double w167 = 0.077751058491018276949*h0*h2;  
     const double w168 = 0.020833333333333333333*h0*h2;  
     const double w169 = 0.0055822748423150563848*h0*h2;  
2216      const double w17 = -0.0009303791403858427308*h0;      const double w17 = -0.0009303791403858427308*h0;
     const double w170 = 0.077751058491018276949*h0*h1;  
     const double w171 = 0.020833333333333333333*h0*h1;  
     const double w172 = 0.0055822748423150563848*h0*h1;  
     const double w173 = -0.25*h1*h2;  
     const double w174 = -0.25*h0*h2;  
     const double w175 = -0.25*h0*h1;  
     const double w176 = 0.25*h1*h2;  
     const double w177 = 0.25*h0*h2;  
     const double w178 = 0.25*h0*h1;  
     const double w179 = 0.061320326520293008568*h0*h1*h2;  
2217      const double w18 = 0.012958509748503046158*h1;      const double w18 = 0.012958509748503046158*h1;
     const double w180 = 0.01643073197072526838*h0*h1*h2;  
     const double w181 = 0.004402601362608064953*h0*h1*h2;  
     const double w182 = 0.0011796734797069914318*h0*h1*h2;  
     const double w183 = 0.125*h0*h1*h2;  
2218      const double w19 = 0.0034722222222222222222*h0;      const double w19 = 0.0034722222222222222222*h0;
     const double w2 = -0.00024929433932114870101*h1;  
2219      const double w20 = 0.012958509748503046158*h2;      const double w20 = 0.012958509748503046158*h2;
2220      const double w21 = -0.012958509748503046158*h1;      const double w21 = -0.012958509748503046158*h1;
2221      const double w22 = -0.012958509748503046158*h0;      const double w22 = -0.012958509748503046158*h0;
# Line 2319  void Brick::assemblePDESingle(Paso_Syste Line 2226  void Brick::assemblePDESingle(Paso_Syste
2226      const double w27 = 0.00024929433932114870101*h0;      const double w27 = 0.00024929433932114870101*h0;
2227      const double w28 = -0.04836181677178996241*h1;      const double w28 = -0.04836181677178996241*h1;
2228      const double w29 = -0.04836181677178996241*h0;      const double w29 = -0.04836181677178996241*h0;
     const double w3 = 0.0009303791403858427308*h0*h2/h1;  
2229      const double w30 = -0.0009303791403858427308*h1*h2/h0;      const double w30 = -0.0009303791403858427308*h1*h2/h0;
2230      const double w31 = -0.0009303791403858427308*h2;      const double w31 = -0.0009303791403858427308*h2;
2231      const double w32 = -0.0009303791403858427308*h0*h2/h1;      const double w32 = -0.0009303791403858427308*h0*h2/h1;
# Line 2330  void Brick::assemblePDESingle(Paso_Syste Line 2236  void Brick::assemblePDESingle(Paso_Syste
2236      const double w37 = -0.012958509748503046158*h2;      const double w37 = -0.012958509748503046158*h2;
2237      const double w38 = -0.012958509748503046158*h0*h2/h1;      const double w38 = -0.012958509748503046158*h0*h2/h1;
2238      const double w39 = -0.04836181677178996241*h2;      const double w39 = -0.04836181677178996241*h2;
     const double w4 = -0.00024929433932114870101*h0;  
2239      const double w40 = -0.0034722222222222222222*h0*h2/h1;      const double w40 = -0.0034722222222222222222*h0*h2/h1;
2240      const double w41 = 0.0009303791403858427308*h0*h1/h2;      const double w41 = 0.0009303791403858427308*h0*h1/h2;
2241      const double w42 = 0.04836181677178996241*h2;      const double w42 = 0.04836181677178996241*h2;
# Line 2341  void Brick::assemblePDESingle(Paso_Syste Line 2246  void Brick::assemblePDESingle(Paso_Syste
2246      const double w47 = -0.0034722222222222222222*h1*h2/h0;      const double w47 = -0.0034722222222222222222*h1*h2/h0;
2247      const double w48 = -0.00024929433932114870101*h1*h2/h0;      const double w48 = -0.00024929433932114870101*h1*h2/h0;
2248      const double w49 = -0.04836181677178996241*h1*h2/h0;      const double w49 = -0.04836181677178996241*h1*h2/h0;
     const double w5 = 0.0009303791403858427308*h1;  
2249      const double w50 = 0.0034722222222222222222*h0*h2/h1;      const double w50 = 0.0034722222222222222222*h0*h2/h1;
2250      const double w51 = -0.0009303791403858427308*h0*h1/h2;      const double w51 = -0.0009303791403858427308*h0*h1/h2;
2251      const double w52 = -0.012958509748503046158*h0*h1/h2;      const double w52 = -0.012958509748503046158*h0*h1/h2;
# Line 2352  void Brick::assemblePDESingle(Paso_Syste Line 2256  void Brick::assemblePDESingle(Paso_Syste
2256      const double w57 = 0.04836181677178996241*h0*h1/h2;      const double w57 = 0.04836181677178996241*h0*h1/h2;
2257      const double w58 = 0.00024929433932114870101*h1*h2/h0;      const double w58 = 0.00024929433932114870101*h1*h2/h0;
2258      const double w59 = 0.00024929433932114870101*h0*h2/h1;      const double w59 = 0.00024929433932114870101*h0*h2/h1;
     const double w6 = 0.0009303791403858427308*h0;  
2259      const double w60 = 0.055555555555555555556*h1*h2/h0;      const double w60 = 0.055555555555555555556*h1*h2/h0;
2260      const double w61 = 0.041666666666666666667*h2;      const double w61 = 0.041666666666666666667*h2;
2261      const double w62 = -0.083333333333333333333*h1;      const double w62 = -0.083333333333333333333*h1;
# Line 2363  void Brick::assemblePDESingle(Paso_Syste Line 2266  void Brick::assemblePDESingle(Paso_Syste
2266      const double w67 = -0.11111111111111111111*h0*h1/h2;      const double w67 = -0.11111111111111111111*h0*h1/h2;
2267      const double w68 = -0.055555555555555555556*h1*h2/h0;      const double w68 = -0.055555555555555555556*h1*h2/h0;
2268      const double w69 = -0.083333333333333333333*h2;      const double w69 = -0.083333333333333333333*h2;
     const double w7 = -0.00024929433932114870101*h0*h1/h2;  
2269      const double w70 = -0.041666666666666666667*h1;      const double w70 = -0.041666666666666666667*h1;
2270      const double w71 = -0.055555555555555555556*h0*h2/h1;      const double w71 = -0.055555555555555555556*h0*h2/h1;
2271      const double w72 = -0.041666666666666666667*h0;      const double w72 = -0.041666666666666666667*h0;
# Line 2374  void Brick::assemblePDESingle(Paso_Syste Line 2276  void Brick::assemblePDESingle(Paso_Syste
2276      const double w77 = -0.11111111111111111111*h0*h2/h1;      const double w77 = -0.11111111111111111111*h0*h2/h1;
2277      const double w78 = 0.055555555555555555556*h0*h1/h2;      const double w78 = 0.055555555555555555556*h0*h1/h2;
2278      const double w79 = -0.11111111111111111111*h1*h2/h0;      const double w79 = -0.11111111111111111111*h1*h2/h0;
     const double w8 = 0.0034722222222222222222*h2;  
2279      const double w80 = -0.027777777777777777778*h1*h2/h0;      const double w80 = -0.027777777777777777778*h1*h2/h0;
2280      const double w81 = -0.041666666666666666667*h2;      const double w81 = -0.041666666666666666667*h2;
2281      const double w82 = -0.027777777777777777778*h0*h2/h1;      const double w82 = -0.027777777777777777778*h0*h2/h1;
# Line 2385  void Brick::assemblePDESingle(Paso_Syste Line 2286  void Brick::assemblePDESingle(Paso_Syste
2286      const double w87 = 0.11111111111111111111*h0*h2/h1;      const double w87 = 0.11111111111111111111*h0*h2/h1;
2287      const double w88 = 0.11111111111111111111*h0*h1/h2;      const double w88 = 0.11111111111111111111*h0*h1/h2;
2288      const double w89 = 0.027777777777777777778*h1*h2/h0;      const double w89 = 0.027777777777777777778*h1*h2/h0;
     const double w9 = -0.0009303791403858427308*h1;  
2289      const double w90 = 0.0001966122466178319053*h1*h2;      const double w90 = 0.0001966122466178319053*h1*h2;
2290      const double w91 = 0.0001966122466178319053*h0*h2;      const double w91 = 0.0001966122466178319053*h0*h2;
2291      const double w92 = 0.0001966122466178319053*h0*h1;      const double w92 = 0.0001966122466178319053*h0*h1;
# Line 2396  void Brick::assemblePDESingle(Paso_Syste Line 2296  void Brick::assemblePDESingle(Paso_Syste
2296      const double w97 = 0.0007337668937680108255*h0*h2;      const double w97 = 0.0007337668937680108255*h0*h2;
2297      const double w98 = 0.010220054420048834761*h1*h2;      const double w98 = 0.010220054420048834761*h1*h2;
2298      const double w99 = 0.010220054420048834761*h0*h2;      const double w99 = 0.010220054420048834761*h0*h2;
2299      /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */      const double w100 = 0.038141762351741127649*h0*h1;
2300        const double w101 = 0.000052682092703316795705*h0*h1;
2301        const double w102 = 0.0007337668937680108255*h0*h1;
2302        const double w103 = 0.010220054420048834761*h0*h1;
2303        const double w104 = -0.0001966122466178319053*h1*h2;
2304        const double w105 = -0.0001966122466178319053*h0*h2;
2305        const double w106 = -0.0007337668937680108255*h1*h2;
2306        const double w107 = -0.0007337668937680108255*h0*h2;
2307        const double w108 = -0.0027384553284542113967*h1*h2;
2308        const double w109 = -0.0027384553284542113967*h0*h2;
2309        const double w110 = -0.010220054420048834761*h1*h2;
2310        const double w111 = -0.010220054420048834761*h0*h2;
2311        const double w112 = -0.0007337668937680108255*h0*h1;
2312        const double w113 = -0.010220054420048834761*h0*h1;
2313        const double w114 = -0.038141762351741127649*h0*h2;
2314        const double w115 = -0.000052682092703316795705*h0*h2;
2315        const double w116 = -0.0001966122466178319053*h0*h1;
2316        const double w117 = -0.0027384553284542113967*h0*h1;
2317        const double w118 = 0.000052682092703316795705*h0*h2;
2318        const double w119 = 0.038141762351741127649*h0*h2;
2319        const double w120 = 0.000052682092703316795705*h1*h2;
2320        const double w121 = 0.038141762351741127649*h1*h2;
2321        const double w122 = -0.000052682092703316795705*h0*h1;
2322        const double w123 = -0.038141762351741127649*h0*h1;
2323        const double w124 = -0.000052682092703316795705*h1*h2;
2324        const double w125 = -0.038141762351741127649*h1*h2;
2325        const double w126 = 0.027777777777777777778*h1*h2;
2326        const double w127 = 0.027777777777777777778*h0*h2;
2327        const double w128 = 0.055555555555555555556*h0*h1;
2328        const double w129 = -0.027777777777777777778*h1*h2;
2329        const double w130 = -0.027777777777777777778*h0*h2;
2330        const double w131 = 0.013888888888888888889*h0*h1;
2331        const double w132 = -0.055555555555555555556*h0*h2;
2332        const double w133 = -0.027777777777777777778*h0*h1;
2333        const double w134 = 0.055555555555555555556*h0*h2;
2334        const double w135 = 0.027777777777777777778*h0*h1;
2335        const double w136 = -0.013888888888888888889*h0*h1;
2336        const double w137 = 0.055555555555555555556*h1*h2;
2337        const double w138 = -0.013888888888888888889*h1*h2;
2338        const double w139 = -0.013888888888888888889*h0*h2;
2339        const double w140 = -0.055555555555555555556*h0*h1;
2340        const double w141 = 0.013888888888888888889*h1*h2;
2341        const double w142 = 0.013888888888888888889*h0*h2;
2342        const double w143 = -0.055555555555555555556*h1*h2;
2343        const double w144 = 0.000041549056553524783501*h0*h1*h2;
2344        const double w145 = 0.0005787037037037037037*h0*h1*h2;
2345        const double w146 = 0.0080603027952983270684*h0*h1*h2;
2346        const double w147 = 0.0001550631900643071218*h0*h1*h2;
2347        const double w148 = 0.002159751624750507693*h0*h1*h2;
2348        const double w149 = 0.03008145955644280058*h0*h1*h2;
2349        const double w150 = 0.000011133036149792012204*h0*h1*h2;
2350        const double w151 = 0.018518518518518518519*h0*h1*h2;
2351        const double w152 = 0.0092592592592592592592*h0*h1*h2;
2352        const double w153 = 0.0046296296296296296296*h0*h1*h2;
2353        const double w154 = 0.037037037037037037037*h0*h1*h2;
2354        const double w155 = -0.077751058491018276949*h1*h2;
2355        const double w156 = -0.077751058491018276949*h0*h2;
2356        const double w157 = -0.077751058491018276949*h0*h1;
2357        const double w158 = -0.020833333333333333333*h0*h2;
2358        const double w159 = -0.020833333333333333333*h0*h1;
2359        const double w160 = -0.020833333333333333333*h1*h2;
2360        const double w161 = -0.0055822748423150563848*h0*h1;
2361        const double w162 = -0.0055822748423150563848*h0*h2;
2362        const double w163 = -0.0055822748423150563848*h1*h2;
2363        const double w164 = 0.077751058491018276949*h1*h2;
2364        const double w165 = 0.020833333333333333333*h1*h2;
2365        const double w166 = 0.0055822748423150563848*h1*h2;
2366        const double w167 = 0.077751058491018276949*h0*h2;
2367        const double w168 = 0.020833333333333333333*h0*h2;
2368        const double w169 = 0.0055822748423150563848*h0*h2;
2369        const double w170 = 0.077751058491018276949*h0*h1;
2370        const double w171 = 0.020833333333333333333*h0*h1;
2371        const double w172 = 0.0055822748423150563848*h0*h1;
2372        const double w173 = -0.25*h1*h2;
2373        const double w174 = -0.25*h0*h2;
2374        const double w175 = -0.25*h0*h1;
2375        const double w176 = 0.25*h1*h2;
2376        const double w177 = 0.25*h0*h2;
2377        const double w178 = 0.25*h0*h1;
2378        const double w179 = 0.061320326520293008568*h0*h1*h2;
2379        const double w180 = 0.01643073197072526838*h0*h1*h2;
2380        const double w181 = 0.004402601362608064953*h0*h1*h2;
2381        const double w182 = 0.0011796734797069914318*h0*h1*h2;
2382        const double w183 = 0.125*h0*h1*h2;
2383    
2384      rhs.requireWrite();      rhs.requireWrite();
2385  #pragma omp parallel  #pragma omp parallel
# Line 2411  void Brick::assemblePDESingle(Paso_Syste Line 2394  void Brick::assemblePDESingle(Paso_Syste
2394                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
2395                          vector<double> EM_F(8, 0);                          vector<double> EM_F(8, 0);
2396                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;
                         /*** GENERATOR SNIP_PDE_SINGLE TOP */  
2397                          ///////////////                          ///////////////
2398                          // process A //                          // process A //
2399                          ///////////////                          ///////////////
# Line 3326  void Brick::assemblePDESingle(Paso_Syste Line 3308  void Brick::assemblePDESingle(Paso_Syste
3308                                  EM_S[INDEX2(0,5,8)]+=tmp118_1 + tmp120_1 + tmp242_1 + tmp248_1 + tmp250_1 + tmp253_1 + tmp334_1 + tmp339_1 + tmp35_1 + tmp37_1 + tmp39_1 + tmp403_1 + tmp410_1 + tmp414_1 + tmp416_1 + tmp417_1 + tmp418_1 + tmp421_1 + tmp422_1 + tmp50_1 + tmp536_1 + tmp537_1;                                  EM_S[INDEX2(0,5,8)]+=tmp118_1 + tmp120_1 + tmp242_1 + tmp248_1 + tmp250_1 + tmp253_1 + tmp334_1 + tmp339_1 + tmp35_1 + tmp37_1 + tmp39_1 + tmp403_1 + tmp410_1 + tmp414_1 + tmp416_1 + tmp417_1 + tmp418_1 + tmp421_1 + tmp422_1 + tmp50_1 + tmp536_1 + tmp537_1;
3309                                  EM_S[INDEX2(3,4,8)]+=tmp155_1 + tmp158_1 + tmp159_1 + tmp161_1 + tmp191_1 + tmp192_1 + tmp193_1 + tmp224_1 + tmp225_1 + tmp226_1 + tmp518_1 + tmp520_1;                                  EM_S[INDEX2(3,4,8)]+=tmp155_1 + tmp158_1 + tmp159_1 + tmp161_1 + tmp191_1 + tmp192_1 + tmp193_1 + tmp224_1 + tmp225_1 + tmp226_1 + tmp518_1 + tmp520_1;
3310                                  EM_S[INDEX2(2,4,8)]+=tmp114_1 + tmp119_1 + tmp240_1 + tmp248_1 + tmp250_1 + tmp254_1 + tmp32_1 + tmp407_1 + tmp411_1 + tmp416_1 + tmp421_1 + tmp42_1 + tmp450_1 + tmp451_1 + tmp47_1 + tmp484_1 + tmp485_1 + tmp487_1 + tmp490_1 + tmp51_1 + tmp538_1 + tmp539_1;                                  EM_S[INDEX2(2,4,8)]+=tmp114_1 + tmp119_1 + tmp240_1 + tmp248_1 + tmp250_1 + tmp254_1 + tmp32_1 + tmp407_1 + tmp411_1 + tmp416_1 + tmp421_1 + tmp42_1 + tmp450_1 + tmp451_1 + tmp47_1 + tmp484_1 + tmp485_1 + tmp487_1 + tmp490_1 + tmp51_1 + tmp538_1 + tmp539_1;
3311                              } else { /* constant data */                              } else { // constant data
3312                                  const register double A_00 = A_p[INDEX2(0,0,3)];                                  const register double A_00 = A_p[INDEX2(0,0,3)];
3313                                  const register double A_01 = A_p[INDEX2(0,1,3)];                                  const register double A_01 = A_p[INDEX2(0,1,3)];
3314                                  const register double A_02 = A_p[INDEX2(0,2,3)];                                  const register double A_02 = A_p[INDEX2(0,2,3)];
# Line 4052  void Brick::assemblePDESingle(Paso_Syste Line 4034  void Brick::assemblePDESingle(Paso_Syste
4034                                  EM_S[INDEX2(0,5,8)]+=tmp290_1 + tmp291_1 + tmp294_1 + tmp295_1 + tmp297_1 + tmp299_1 + tmp384_1 + tmp385_1 + tmp386_1 + tmp387_1;                                  EM_S[INDEX2(0,5,8)]+=tmp290_1 + tmp291_1 + tmp294_1 + tmp295_1 + tmp297_1 + tmp299_1 + tmp384_1 + tmp385_1 + tmp386_1 + tmp387_1;
4035                                  EM_S[INDEX2(3,4,8)]+=tmp100_1 + tmp101_1 + tmp122_1 + tmp123_1 + tmp80_1 + tmp81_1;                                  EM_S[INDEX2(3,4,8)]+=tmp100_1 + tmp101_1 + tmp122_1 + tmp123_1 + tmp80_1 + tmp81_1;
4036                                  EM_S[INDEX2(2,4,8)]+=tmp388_1 + tmp389_1 + tmp392_1 + tmp394_1 + tmp395_1 + tmp396_1 + tmp448_1 + tmp449_1 + tmp450_1 + tmp451_1;                                  EM_S[INDEX2(2,4,8)]+=tmp388_1 + tmp389_1 + tmp392_1 + tmp394_1 + tmp395_1 + tmp396_1 + tmp448_1 + tmp449_1 + tmp450_1 + tmp451_1;
4037                              } else { /* constant data */                              } else { // constant data
4038                                  const register double B_0 = B_p[0];                                  const register double B_0 = B_p[0];
4039                                  const register double B_1 = B_p[1];                                  const register double B_1 = B_p[1];
4040                                  const register double B_2 = B_p[2];                                  const register double B_2 = B_p[2];
# Line 4733  void Brick::assemblePDESingle(Paso_Syste Line 4715  void Brick::assemblePDESingle(Paso_Syste
4715                                  EM_S[INDEX2(0,5,8)]+=tmp290_1 + tmp291_1 + tmp292_1 + tmp293_1 + tmp296_1 + tmp299_1 + tmp384_1 + tmp385_1 + tmp386_1 + tmp387_1;                                  EM_S[INDEX2(0,5,8)]+=tmp290_1 + tmp291_1 + tmp292_1 + tmp293_1 + tmp296_1 + tmp299_1 + tmp384_1 + tmp385_1 + tmp386_1 + tmp387_1;
4716                                  EM_S[INDEX2(3,4,8)]+=tmp100_1 + tmp101_1 + tmp122_1 + tmp123_1 + tmp78_1 + tmp81_1;                                  EM_S[INDEX2(3,4,8)]+=tmp100_1 + tmp101_1 + tmp122_1 + tmp123_1 + tmp78_1 + tmp81_1;
4717                                  EM_S[INDEX2(2,4,8)]+=tmp388_1 + tmp389_1 + tmp390_1 + tmp391_1 + tmp396_1 + tmp397_1 + tmp448_1 + tmp449_1 + tmp450_1 + tmp451_1;                                  EM_S[INDEX2(2,4,8)]+=tmp388_1 + tmp389_1 + tmp390_1 + tmp391_1 + tmp396_1 + tmp397_1 + tmp448_1 + tmp449_1 + tmp450_1 + tmp451_1;
4718                              } else { /* constant data */                              } else { // constant data
4719                                  const register double C_0 = C_p[0];                                  const register double C_0 = C_p[0];
4720                                  const register double C_1 = C_p[1];                                  const register double C_1 = C_p[1];
4721                                  const register double C_2 = C_p[2];                                  const register double C_2 = C_p[2];
# Line 5008  void Brick::assemblePDESingle(Paso_Syste Line 4990  void Brick::assemblePDESingle(Paso_Syste
4990                                  EM_S[INDEX2(5,7,8)]+=tmp32_1 + tmp47_1 + tmp48_1;                                  EM_S[INDEX2(5,7,8)]+=tmp32_1 + tmp47_1 + tmp48_1;
4991                                  EM_S[INDEX2(6,7,8)]+=tmp25_1 + tmp26_1 + tmp27_1;                                  EM_S[INDEX2(6,7,8)]+=tmp25_1 + tmp26_1 + tmp27_1;
4992                                  EM_S[INDEX2(7,7,8)]+=tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1;                                  EM_S[INDEX2(7,7,8)]+=tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1;
4993                              } else { /* constant data */                              } else { // constant data
4994                                  const register double D_0 = D_p[0];                                  const register double D_0 = D_p[0];
4995                                  const register double tmp3_1 = D_0*w154;                                  const register double tmp3_1 = D_0*w154;
4996                                  const register double tmp0_1 = D_0*w151;                                  const register double tmp0_1 = D_0*w151;
# Line 5197  void Brick::assemblePDESingle(Paso_Syste Line 5179  void Brick::assemblePDESingle(Paso_Syste
5179                                  EM_F[5]+=tmp29_1 + tmp41_1 + tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1 + tmp7_1;                                  EM_F[5]+=tmp29_1 + tmp41_1 + tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1 + tmp7_1;
5180                                  EM_F[6]+=tmp28_1 + tmp2_1 + tmp43_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1 + tmp52_1 + tmp53_1;                                  EM_F[6]+=tmp28_1 + tmp2_1 + tmp43_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1 + tmp52_1 + tmp53_1;
5181                                  EM_F[7]+=tmp12_1 + tmp21_1 + tmp37_1 + tmp54_1 + tmp55_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;                                  EM_F[7]+=tmp12_1 + tmp21_1 + tmp37_1 + tmp54_1 + tmp55_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
5182                              } else { /* constant data */                              } else { // constant data
5183                                  const register double X_0 = X_p[0];                                  const register double X_0 = X_p[0];
5184                                  const register double X_1 = X_p[1];                                  const register double X_1 = X_p[1];
5185                                  const register double X_2 = X_p[2];                                  const register double X_2 = X_p[2];
# Line 5280  void Brick::assemblePDESingle(Paso_Syste Line 5262  void Brick::assemblePDESingle(Paso_Syste
5262                                  EM_F[5]+=tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;                                  EM_F[5]+=tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;
5263                                  EM_F[6]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1;                                  EM_F[6]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1;
5264                                  EM_F[7]+=tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                                  EM_F[7]+=tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
5265                              } else { /* constant data */                              } else { // constant data
5266                                  const register double Y_0 = Y_p[0];                                  const register double Y_0 = Y_p[0];
5267                                  const register double tmp0_1 = Y_0*w183;                                  const register double tmp0_1 = Y_0*w183;
5268                                  EM_F[0]+=tmp0_1;                                  EM_F[0]+=tmp0_1;
# Line 5293  void Brick::assemblePDESingle(Paso_Syste Line 5275  void Brick::assemblePDESingle(Paso_Syste
5275                                  EM_F[7]+=tmp0_1;                                  EM_F[7]+=tmp0_1;
5276                              }                              }
5277                          }                          }
                         /* GENERATOR SNIP_PDE_SINGLE BOTTOM */  
5278    
5279                          // 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)
5280                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
# Line 5307  void Brick::assemblePDESingle(Paso_Syste Line 5288  void Brick::assemblePDESingle(Paso_Syste
5288                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);
5289                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);
5290                          if (add_EM_F) {                          if (add_EM_F) {
                             //cout << "-- AddtoRHS -- " << endl;  
5291                              double *F_p=rhs.getSampleDataRW(0);                              double *F_p=rhs.getSampleDataRW(0);
5292                              for (index_t i=0; i<rowIndex.size(); i++) {                              for (index_t i=0; i<rowIndex.size(); i++) {
5293                                  if (rowIndex[i]<getNumDOF()) {                                  if (rowIndex[i]<getNumDOF()) {
5294                                      F_p[rowIndex[i]]+=EM_F[i];                                      F_p[rowIndex[i]]+=EM_F[i];
                                     //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;  
5295                                  }                                  }
5296                              }                              }
                             //cout << "---"<<endl;  
5297                          }                          }
5298                          if (add_EM_S) {                          if (add_EM_S) {
                             //cout << "-- AddtoSystem -- " << endl;  
5299                              addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);                              addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
5300                          }                          }
5301                      } // end k0 loop                      } // end k0 loop
# Line 5338  void Brick::assemblePDESingleReduced(Pas Line 5315  void Brick::assemblePDESingleReduced(Pas
5315      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
5316      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
5317      const double h2 = m_l2/m_gNE2;      const double h2 = m_l2/m_gNE2;
     /* GENERATOR SNIP_PDE_SINGLE_REDUCED_PRE TOP */  
5318      const double w0 = 0.0625*h1*h2/h0;      const double w0 = 0.0625*h1*h2/h0;
5319      const double w1 = 0.0625*h2;      const double w1 = 0.0625*h2;
5320      const double w2 = -0.0625*h1;      const double w2 = -0.0625*h1;
# Line 5365  void Brick::assemblePDESingleReduced(Pas Line 5341  void Brick::assemblePDESingleReduced(Pas
5341      const double w23 = 0.25*h0*h2;      const double w23 = 0.25*h0*h2;
5342      const double w24 = 0.25*h0*h1;      const double w24 = 0.25*h0*h1;
5343      const double w25 = 0.125*h0*h1*h2;      const double w25 = 0.125*h0*h1*h2;
     /* GENERATOR SNIP_PDE_SINGLE_REDUCED_PRE BOTTOM */  
5344    
5345      rhs.requireWrite();      rhs.requireWrite();
5346  #pragma omp parallel  #pragma omp parallel
# Line 5380  void Brick::assemblePDESingleReduced(Pas Line 5355  void Brick::assemblePDESingleReduced(Pas
5355                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
5356                          vector<double> EM_F(8, 0);                          vector<double> EM_F(8, 0);
5357                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;
                         /* GENERATOR SNIP_PDE_SINGLE_REDUCED TOP */  
5358                          ///////////////                          ///////////////
5359                          // process A //                          // process A //
5360                          ///////////////                          ///////////////
# Line 5730  void Brick::assemblePDESingleReduced(Pas Line 5704  void Brick::assemblePDESingleReduced(Pas
5704                              const register double X_0 = X_p[0];                              const register double X_0 = X_p[0];
5705                              const register double X_1 = X_p[1];                              const register double X_1 = X_p[1];
5706                              const register double X_2 = X_p[2];                              const register double X_2 = X_p[2];
5707                                const register double tmp0_1 = X_2*w21;
5708                              const register double tmp1_1 = X_0*w19;                              const register double tmp1_1 = X_0*w19;
5709                              const register double tmp2_1 = X_1*w20;                              const register double tmp2_1 = X_1*w20;
5710                              const register double tmp3_1 = X_0*w22;                              const register double tmp3_1 = X_0*w22;
5711                              const register double tmp4_1 = X_1*w23;                              const register double tmp4_1 = X_1*w23;
5712                              const register double tmp5_1 = X_2*w24;                              const register double tmp5_1 = X_2*w24;
                             const register double tmp0_1 = X_2*w21;  
5713                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;                              EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
5714                              EM_F[1]+=tmp0_1 + tmp2_1 + tmp3_1;                              EM_F[1]+=tmp0_1 + tmp2_1 + tmp3_1;
5715                              EM_F[2]+=tmp0_1 + tmp1_1 + tmp4_1;                              EM_F[2]+=tmp0_1 + tmp1_1 + tmp4_1;
# Line 5751  void Brick::assemblePDESingleReduced(Pas Line 5725  void Brick::assemblePDESingleReduced(Pas
5725                          if (!Y.isEmpty()) {                          if (!Y.isEmpty()) {
5726                              add_EM_F=true;                              add_EM_F=true;
5727                              const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                              const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
5728                              const register double Y_0 = Y_p[0];                              const register double tmp0_1 = Y_p[0]*w25;
                             const register double tmp0_1 = Y_0*w25;  
5729                              EM_F[0]+=tmp0_1;                              EM_F[0]+=tmp0_1;
5730                              EM_F[1]+=tmp0_1;                              EM_F[1]+=tmp0_1;
5731                              EM_F[2]+=tmp0_1;                              EM_F[2]+=tmp0_1;
# Line 5762  void Brick::assemblePDESingleReduced(Pas Line 5735  void Brick::assemblePDESingleReduced(Pas
5735                              EM_F[6]+=tmp0_1;                              EM_F[6]+=tmp0_1;
5736                              EM_F[7]+=tmp0_1;                              EM_F[7]+=tmp0_1;
5737                          }                          }
                         /* GENERATOR SNIP_PDE_SINGLE_REDUCED BOTTOM */  
5738    
5739                          // 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)
5740                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
# Line 5776  void Brick::assemblePDESingleReduced(Pas Line 5748  void Brick::assemblePDESingleReduced(Pas
5748                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);
5749                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);
5750                          if (add_EM_F) {                          if (add_EM_F) {
                             //cout << "-- AddtoRHS -- " << endl;  
5751                              double *F_p=rhs.getSampleDataRW(0);                              double *F_p=rhs.getSampleDataRW(0);
5752                              for (index_t i=0; i<rowIndex.size(); i++) {                              for (index_t i=0; i<rowIndex.size(); i++) {
5753                                  if (rowIndex[i]<getNumDOF()) {                                  if (rowIndex[i]<getNumDOF()) {
5754                                      F_p[rowIndex[i]]+=EM_F[i];                                      F_p[rowIndex[i]]+=EM_F[i];
                                     //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;  
5755                                  }                                  }
5756                              }                              }
                             //cout << "---"<<endl;  
5757                          }                          }
5758                          if (add_EM_S) {                          if (add_EM_S) {
                             //cout << "-- AddtoSystem -- " << endl;  
5759                              addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);                              addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
5760                          }                          }
5761                      } // end k0 loop                      } // end k0 loop
# Line 5815  void Brick::assemblePDESystem(Paso_Syste Line 5783  void Brick::assemblePDESystem(Paso_Syste
5783          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
5784      }      }
5785    
     /*** GENERATOR SNIP_PDE_SYSTEM_PRE TOP */  
5786      const double w0 = 0.0009303791403858427308*h1*h2/h0;      const double w0 = 0.0009303791403858427308*h1*h2/h0;
5787      const double w1 = 0.0009303791403858427308*h2;      const double w1 = 0.0009303791403858427308*h2;
5788        const double w2 = -0.00024929433932114870101*h1;
5789        const double w3 = 0.0009303791403858427308*h0*h2/h1;
5790        const double w4 = -0.00024929433932114870101*h0;
5791        const double w5 = 0.0009303791403858427308*h1;
5792        const double w6 = 0.0009303791403858427308*h0;
5793        const double w7 = -0.00024929433932114870101*h0*h1/h2;
5794        const double w8 = 0.0034722222222222222222*h2;
5795        const double w9 = -0.0009303791403858427308*h1;
5796      const double w10 = 0.012958509748503046158*h0*h2/h1;      const double w10 = 0.012958509748503046158*h0*h2/h1;
     const double w100 = 0.038141762351741127649*h0*h1;  
     const double w101 = 0.000052682092703316795705*h0*h1;  
     const double w102 = 0.0007337668937680108255*h0*h1;  
     const double w103 = 0.010220054420048834761*h0*h1;  
     const double w104 = -0.0001966122466178319053*h1*h2;  
     const double w105 = -0.0001966122466178319053*h0*h2;  
     const double w106 = -0.0007337668937680108255*h1*h2;  
     const double w107 = -0.0007337668937680108255*h0*h2;  
     const double w108 = -0.0027384553284542113967*h1*h2;  
     const double w109 = -0.0027384553284542113967*h0*h2;  
5797      const double w11 = -0.0034722222222222222222*h0;      const double w11 = -0.0034722222222222222222*h0;
     const double w110 = -0.010220054420048834761*h1*h2;  
     const double w111 = -0.010220054420048834761*h0*h2;  
     const double w112 = -0.0007337668937680108255*h0*h1;  
     const double w113 = -0.010220054420048834761*h0*h1;  
     const double w114 = -0.038141762351741127649*h0*h2;  
     const double w115 = -0.000052682092703316795705*h0*h2;  
     const double w116 = -0.0001966122466178319053*h0*h1;  
     const double w117 = -0.0027384553284542113967*h0*h1;  
     const double w118 = 0.000052682092703316795705*h0*h2;  
     const double w119 = 0.038141762351741127649*h0*h2;  
5798      const double w12 = 0.0034722222222222222222*h1;      const double w12 = 0.0034722222222222222222*h1;
     const double w120 = 0.000052682092703316795705*h1*h2;  
     const double w121 = 0.038141762351741127649*h1*h2;  
     const double w122 = -0.000052682092703316795705*h0*h1;  
     const double w123 = -0.038141762351741127649*h0*h1;  
     const double w124 = -0.000052682092703316795705*h1*h2;  
     const double w125 = -0.038141762351741127649*h1*h2;  
     const double w126 = 0.027777777777777777778*h1*h2;  
     const double w127 = 0.027777777777777777778*h0*h2;  
     const double w128 = 0.055555555555555555556*h0*h1;  
     const double w129 = -0.027777777777777777778*h1*h2;  
5799      const double w13 = 0.012958509748503046158*h0;      const double w13 = 0.012958509748503046158*h0;
     const double w130 = -0.027777777777777777778*h0*h2;  
     const double w131 = 0.013888888888888888889*h0*h1;  
     const double w132 = -0.055555555555555555556*h0*h2;  
     const double w133 = -0.027777777777777777778*h0*h1;  
     const double w134 = 0.055555555555555555556*h0*h2;  
     const double w135 = 0.027777777777777777778*h0*h1;  
     const double w136 = -0.013888888888888888889*h0*h1;  
     const double w137 = 0.055555555555555555556*h1*h2;  
     const double w138 = -0.013888888888888888889*h1*h2;  
     const double w139 = -0.013888888888888888889*h0*h2;  
5800      const double w14 = -0.0034722222222222222222*h0*h1/h2;      const double w14 = -0.0034722222222222222222*h0*h1/h2;
     const double w140 = -0.055555555555555555556*h0*h1;  
     const double w141 = 0.013888888888888888889*h1*h2;  
     const double w142 = 0.013888888888888888889*h0*h2;  
     const double w143 = -0.055555555555555555556*h1*h2;  
     const double w144 = 0.000041549056553524783501*h0*h1*h2;  
     const double w145 = 0.0005787037037037037037*h0*h1*h2;  
     const double w146 = 0.0080603027952983270684*h0*h1*h2;  
     const double w147 = 0.0001550631900643071218*h0*h1*h2;  
     const double w148 = 0.002159751624750507693*h0*h1*h2;  
     const double w149 = 0.03008145955644280058*h0*h1*h2;  
5801      const double w15 = 0.012958509748503046158*h1*h2/h0;      const double w15 = 0.012958509748503046158*h1*h2/h0;
     const double w150 = 0.000011133036149792012204*h0*h1*h2;  
     const double w151 = 0.018518518518518518519*h0*h1*h2;  
     const double w152 = 0.0092592592592592592592*h0*h1*h2;  
     const double w153 = 0.0046296296296296296296*h0*h1*h2;  
     const double w154 = 0.037037037037037037037*h0*h1*h2;  
     const double w155 = -0.077751058491018276949*h1*h2;  
     const double w156 = -0.077751058491018276949*h0*h2;  
     const double w157 = -0.077751058491018276949*h0*h1;  
     const double w158 = -0.020833333333333333333*h0*h2;  
     const double w159 = -0.020833333333333333333*h0*h1;  
5802      const double w16 = -0.0034722222222222222222*h1;      const double w16 = -0.0034722222222222222222*h1;
     const double w160 = -0.020833333333333333333*h1*h2;  
     const double w161 = -0.0055822748423150563848*h0*h1;  
     const double w162 = -0.0055822748423150563848*h0*h2;  
     const double w163 = -0.0055822748423150563848*h1*h2;  
     const double w164 = 0.077751058491018276949*h1*h2;  
     const double w165 = 0.020833333333333333333*h1*h2;  
     const double w166 = 0.0055822748423150563848*h1*h2;  
     const double w167 = 0.077751058491018276949*h0*h2;  
     const double w168 = 0.020833333333333333333*h0*h2;  
     const double w169 = 0.0055822748423150563848*h0*h2;  
5803      const double w17 = -0.0009303791403858427308*h0;      const double w17 = -0.0009303791403858427308*h0;
     const double w170 = 0.077751058491018276949*h0*h1;  
     const double w171 = 0.020833333333333333333*h0*h1;  
     const double w172 = 0.0055822748423150563848*h0*h1;  
     const double w173 = -0.25*h1*h2;  
     const double w174 = -0.25*h0*h2;  
     const double w175 = -0.25*h0*h1;  
     const double w176 = 0.25*h1*h2;  
     const double w177 = 0.25*h0*h2;  
     const double w178 = 0.25*h0*h1;  
     const double w179 = 0.061320326520293008568*h0*h1*h2;  
5804      const double w18 = 0.012958509748503046158*h1;      const double w18 = 0.012958509748503046158*h1;
     const double w180 = 0.01643073197072526838*h0*h1*h2;  
     const double w181 = 0.004402601362608064953*h0*h1*h2;  
     const double w182 = 0.0011796734797069914318*h0*h1*h2;  
     const double w183 = 0.125*h0*h1*h2;  
5805      const double w19 = 0.0034722222222222222222*h0;      const double w19 = 0.0034722222222222222222*h0;
     const double w2 = -0.00024929433932114870101*h1;  
5806      const double w20 = 0.012958509748503046158*h2;      const double w20 = 0.012958509748503046158*h2;
5807      const double w21 = -0.012958509748503046158*h1;      const double w21 = -0.012958509748503046158*h1;
5808      const double w22 = -0.012958509748503046158*h0;      const double w22 = -0.012958509748503046158*h0;
# Line 5923  void Brick::assemblePDESystem(Paso_Syste Line 5813  void Brick::assemblePDESystem(Paso_Syste
5813      const double w27 = 0.00024929433932114870101*h0;      const double w27 = 0.00024929433932114870101*h0;
5814      const double w28 = -0.04836181677178996241*h1;      const double w28 = -0.04836181677178996241*h1;
5815      const double w29 = -0.04836181677178996241*h0;      const double w29 = -0.04836181677178996241*h0;
     const double w3 = 0.0009303791403858427308*h0*h2/h1;  
5816      const double w30 = -0.0009303791403858427308*h1*h2/h0;      const double w30 = -0.0009303791403858427308*h1*h2/h0;
5817      const double w31 = -0.0009303791403858427308*h2;      const double w31 = -0.0009303791403858427308*h2;
5818      const double w32 = -0.0009303791403858427308*h0*h2/h1;      const double w32 = -0.0009303791403858427308*h0*h2/h1;
# Line 5934  void Brick::assemblePDESystem(Paso_Syste Line 5823  void Brick::assemblePDESystem(Paso_Syste
5823      const double w37 = -0.012958509748503046158*h2;      const double w37 = -0.012958509748503046158*h2;
5824      const double w38 = -0.012958509748503046158*h0*h2/h1;      const double w38 = -0.012958509748503046158*h0*h2/h1;
5825      const double w39 = -0.04836181677178996241*h2;      const double w39 = -0.04836181677178996241*h2;
     const double w4 = -0.00024929433932114870101*h0;  
5826      const double w40 = -0.0034722222222222222222*h0*h2/h1;      const double w40 = -0.0034722222222222222222*h0*h2/h1;
5827      const double w41 = 0.0009303791403858427308*h0*h1/h2;      const double w41 = 0.0009303791403858427308*h0*h1/h2;
5828      const double w42 = 0.04836181677178996241*h2;      const double w42 = 0.04836181677178996241*h2;
# Line 5945  void Brick::assemblePDESystem(Paso_Syste Line 5833  void Brick::assemblePDESystem(Paso_Syste
5833      const double w47 = -0.0034722222222222222222*h1*h2/h0;      const double w47 = -0.0034722222222222222222*h1*h2/h0;
5834      const double w48 = -0.00024929433932114870101*h1*h2/h0;      const double w48 = -0.00024929433932114870101*h1*h2/h0;
5835      const double w49 = -0.04836181677178996241*h1*h2/h0;      const double w49 = -0.04836181677178996241*h1*h2/h0;
     const double w5 = 0.0009303791403858427308*h1;  
5836      const double w50 = 0.0034722222222222222222*h0*h2/h1;      const double w50 = 0.0034722222222222222222*h0*h2/h1;
5837      const double w51 = -0.0009303791403858427308*h0*h1/h2;      const double w51 = -0.0009303791403858427308*h0*h1/h2;
5838      const double w52 = -0.012958509748503046158*h0*h1/h2;      const double w52 = -0.012958509748503046158*h0*h1/h2;
# Line 5956  void Brick::assemblePDESystem(Paso_Syste Line 5843  void Brick::assemblePDESystem(Paso_Syste
5843      const double w57 = 0.04836181677178996241*h0*h1/h2;      const double w57 = 0.04836181677178996241*h0*h1/h2;
5844      const double w58 = 0.00024929433932114870101*h1*h2/h0;      const double w58 = 0.00024929433932114870101*h1*h2/h0;
5845      const double w59 = 0.00024929433932114870101*h0*h2/h1;      const double w59 = 0.00024929433932114870101*h0*h2/h1;
     const double w6 = 0.0009303791403858427308*h0;  
5846      const double w60 = 0.055555555555555555556*h1*h2/h0;      const double w60 = 0.055555555555555555556*h1*h2/h0;
5847      const double w61 = 0.041666666666666666667*h2;      const double w61 = 0.041666666666666666667*h2;
5848      const double w62 = -0.083333333333333333333*h1;      const double w62 = -0.083333333333333333333*h1;
# Line 5967  void Brick::assemblePDESystem(Paso_Syste Line 5853  void Brick::assemblePDESystem(Paso_Syste
5853      const double w67 = -0.11111111111111111111*h0*h1/h2;      const double w67 = -0.11111111111111111111*h0*h1/h2;
5854      const double w68 = -0.055555555555555555556*h1*h2/h0;      const double w68 = -0.055555555555555555556*h1*h2/h0;
5855      const double w69 = -0.083333333333333333333*h2;      const double w69 = -0.083333333333333333333*h2;
     const double w7 = -0.00024929433932114870101*h0*h1/h2;  
5856      const double w70 = -0.041666666666666666667*h1;      const double w70 = -0.041666666666666666667*h1;
5857      const double w71 = -0.055555555555555555556*h0*h2/h1;      const double w71 = -0.055555555555555555556*h0*h2/h1;
5858      const double w72 = -0.041666666666666666667*h0;      const double w72 = -0.041666666666666666667*h0;
# Line 5978  void Brick::assemblePDESystem(Paso_Syste Line 5863  void Brick::assemblePDESystem(Paso_Syste
5863      const double w77 = -0.11111111111111111111*h0*h2/h1;      const double w77 = -0.11111111111111111111*h0*h2/h1;
5864      const double w78 = 0.055555555555555555556*h0*h1/h2;      const double w78 = 0.055555555555555555556*h0*h1/h2;
5865      const double w79 = -0.11111111111111111111*h1*h2/h0;      const double w79 = -0.11111111111111111111*h1*h2/h0;
     const double w8 = 0.0034722222222222222222*h2;  
5866      const double w80 = -0.027777777777777777778*h1*h2/h0;      const double w80 = -0.027777777777777777778*h1*h2/h0;
5867      const double w81 = -0.041666666666666666667*h2;      const double w81 = -0.041666666666666666667*h2;
5868      const double w82 = -0.027777777777777777778*h0*h2/h1;      const double w82 = -0.027777777777777777778*h0*h2/h1;
# Line 5989  void Brick::assemblePDESystem(Paso_Syste Line 5873  void Brick::assemblePDESystem(Paso_Syste
5873      const double w87 = 0.11111111111111111111*h0*h2/h1;      const double w87 = 0.11111111111111111111*h0*h2/h1;
5874      const double w88 = 0.11111111111111111111*h0*h1/h2;      const double w88 = 0.11111111111111111111*h0*h1/h2;
5875      const double w89 = 0.027777777777777777778*h1*h2/h0;      const double w89 = 0.027777777777777777778*h1*h2/h0;
     const double w9 = -0.0009303791403858427308*h1;  
5876      const double w90 = 0.0001966122466178319053*h1*h2;      const double w90 = 0.0001966122466178319053*h1*h2;
5877      const double w91 = 0.0001966122466178319053*h0*h2;      const double w91 = 0.0001966122466178319053*h0*h2;
5878      const double w92 = 0.0001966122466178319053*h0*h1;      const double w92 = 0.0001966122466178319053*h0*h1;
# Line 6000  void Brick::assemblePDESystem(Paso_Syste Line 5883  void Brick::assemblePDESystem(Paso_Syste
5883      const double w97 = 0.0007337668937680108255*h0*h2;      const double w97 = 0.0007337668937680108255*h0*h2;
5884      const double w98 = 0.010220054420048834761*h1*h2;      const double w98 = 0.010220054420048834761*h1*h2;
5885      const double w99 = 0.010220054420048834761*h0*h2;      const double w99 = 0.010220054420048834761*h0*h2;
5886      /* GENERATOR SNIP_PDE_SYSTEM_PRE BOTTOM */      const double w100 = 0.038141762351741127649*h0*h1;
5887        const double w101 = 0.000052682092703316795705*h0*h1;
5888        const double w102 = 0.0007337668937680108255*h0*h1;
5889        const double w103 = 0.010220054420048834761*h0*h1;
5890        const double w104 = -0.0001966122466178319053*h1*h2;
5891        const double w105 = -0.0001966122466178319053*h0*h2;
5892        const double w106 = -0.0007337668937680108255*h1*h2;
5893        const double w107 = -0.0007337668937680108255*h0*h2;
5894        const double w108 = -0.0027384553284542113967*h1*h2;
5895        const double w109 = -0.0027384553284542113967*h0*h2;
5896        const double w110 = -0.010220054420048834761*h1*h2;
5897        const double w111 = -0.010220054420048834761*h0*h2;
5898        const double w112 = -0.0007337668937680108255*h0*h1;
5899        const double w113 = -0.010220054420048834761*h0*h1;
5900        const double w114 = -0.038141762351741127649*h0*h2;
5901        const double w115 = -0.000052682092703316795705*h0*h2;
5902        const double w116 = -0.0001966122466178319053*h0*h1;
5903        const double w117 = -0.0027384553284542113967*h0*h1;
5904        const double w118 = 0.000052682092703316795705*h0*h2;
5905        const double w119 = 0.038141762351741127649*h0*h2;
5906        const double w120 = 0.000052682092703316795705*h1*h2;
5907        const double w121 = 0.038141762351741127649*h1*h2;
5908        const double w122 = -0.000052682092703316795705*h0*h1;
5909        const double w123 = -0.038141762351741127649*h0*h1;
5910        const double w124 = -0.000052682092703316795705*h1*h2;
5911        const double w125 = -0.038141762351741127649*h1*h2;
5912        const double w126 = 0.027777777777777777778*h1*h2;
5913        const double w127 = 0.027777777777777777778*h0*h2;
5914        const double w128 = 0.055555555555555555556*h0*h1;
5915        const double w129 = -0.027777777777777777778*h1*h2;
5916        const double w130 = -0.027777777777777777778*h0*h2;
5917        const double w131 = 0.013888888888888888889*h0*h1;
5918        const double w132 = -0.055555555555555555556*h0*h2;
5919        const double w133 = -0.027777777777777777778*h0*h1;
5920        const double w134 = 0.055555555555555555556*h0*h2;
5921        const double w135 = 0.027777777777777777778*h0*h1;
5922        const double w136 = -0.013888888888888888889*h0*h1;
5923        const double w137 = 0.055555555555555555556*h1*h2;
5924        const double w138 = -0.013888888888888888889*h1*h2;
5925        const double w139 = -0.013888888888888888889*h0*h2;
5926        const double w140 = -0.055555555555555555556*h0*h1;
5927        const double w141 = 0.013888888888888888889*h1*h2;
5928        const double w142 = 0.013888888888888888889*h0*h2;
5929        const double w143 = -0.055555555555555555556*h1*h2;
5930        const double w144 = 0.000041549056553524783501*h0*h1*h2;
5931        const double w145 = 0.0005787037037037037037*h0*h1*h2;
5932        const double w146 = 0.0080603027952983270684*h0*h1*h2;
5933        const double w147 = 0.0001550631900643071218*h0*h1*h2;
5934        const double w148 = 0.002159751624750507693*h0*h1*h2;
5935        const double w149 = 0.03008145955644280058*h0*h1*h2;
5936        const double w150 = 0.000011133036149792012204*h0*h1*h2;
5937        const double w151 = 0.018518518518518518519*h0*h1*h2;
5938        const double w152 = 0.0092592592592592592592*h0*h1*h2;
5939        const double w153 = 0.0046296296296296296296*h0*h1*h2;
5940        const double w154 = 0.037037037037037037037*h0*h1*h2;
5941        const double w155 = -0.077751058491018276949*h1*h2;
5942        const double w156 = -0.077751058491018276949*h0*h2;
5943        const double w157 = -0.077751058491018276949*h0*h1;
5944        const double w158 = -0.020833333333333333333*h0*h2;
5945        const double w159 = -0.020833333333333333333*h0*h1;
5946        const double w160 = -0.020833333333333333333*h1*h2;
5947        const double w161 = -0.0055822748423150563848*h0*h1;
5948        const double w162 = -0.0055822748423150563848*h0*h2;
5949        const double w163 = -0.0055822748423150563848*h1*h2;
5950        const double w164 = 0.077751058491018276949*h1*h2;
5951        const double w165 = 0.020833333333333333333*h1*h2;
5952        const double w166 = 0.0055822748423150563848*h1*h2;
5953        const double w167 = 0.077751058491018276949*h0*h2;
5954        const double w168 = 0.020833333333333333333*h0*h2;
5955        const double w169 = 0.0055822748423150563848*h0*h2;
5956        const double w170 = 0.077751058491018276949*h0*h1;
5957        const double w171 = 0.020833333333333333333*h0*h1;
5958        const double w172 = 0.0055822748423150563848*h0*h1;
5959        const double w173 = -0.25*h1*h2;
5960        const double w174 = -0.25*h0*h2;
5961        const double w175 = -0.25*h0*h1;
5962        const double w176 = 0.25*h1*h2;
5963        const double w177 = 0.25*h0*h2;
5964        const double w178 = 0.25*h0*h1;
5965        const double w179 = 0.061320326520293008568*h0*h1*h2;
5966        const double w180 = 0.01643073197072526838*h0*h1*h2;
5967        const double w181 = 0.004402601362608064953*h0*h1*h2;
5968        const double w182 = 0.0011796734797069914318*h0*h1*h2;
5969        const double w183 = 0.125*h0*h1*h2;
5970    
5971      rhs.requireWrite();      rhs.requireWrite();
5972  #pragma omp parallel  #pragma omp parallel
# Line 6015  void Brick::assemblePDESystem(Paso_Syste Line 5981  void Brick::assemblePDESystem(Paso_Syste
5981                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
5982                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
5983                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;
                         /*** GENERATOR SNIP_PDE_SYSTEM TOP */  
5984                          ///////////////                          ///////////////
5985                          // process A //                          // process A //
5986                          ///////////////                          ///////////////
# Line 6934  void Brick::assemblePDESystem(Paso_Syste Line 6899  void Brick::assemblePDESystem(Paso_Syste
6899                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp114_1 + tmp119_1 + tmp240_1 + tmp248_1 + tmp250_1 + tmp254_1 + tmp32_1 + tmp407_1 + tmp411_1 + tmp416_1 + tmp421_1 + tmp42_1 + tmp450_1 + tmp451_1 + tmp47_1 + tmp484_1 + tmp485_1 + tmp487_1 + tmp490_1 + tmp51_1 + tmp538_1 + tmp539_1;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp114_1 + tmp119_1 + tmp240_1 + tmp248_1 + tmp250_1 + tmp254_1 + tmp32_1 + tmp407_1 + tmp411_1 + tmp416_1 + tmp421_1 + tmp42_1 + tmp450_1 + tmp451_1 + tmp47_1 + tmp484_1 + tmp485_1 + tmp487_1 + tmp490_1 + tmp51_1 + tmp538_1 + tmp539_1;
6900                                      }                                      }
6901                                  }                                  }
6902                              } else { /* constant data */                              } else { // constant data
6903                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
6904                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
6905                                          const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,3, numComp)];                                          const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,3, numComp)];
# Line 7668  void Brick::assemblePDESystem(Paso_Syste Line 7633  void Brick::assemblePDESystem(Paso_Syste
7633                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp388_1 + tmp389_1 + tmp392_1 + tmp394_1 + tmp395_1 + tmp396_1 + tmp448_1 + tmp449_1 + tmp450_1 + tmp451_1;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp388_1 + tmp389_1 + tmp392_1 + tmp394_1 + tmp395_1 + tmp396_1 + tmp448_1 + tmp449_1 + tmp450_1 + tmp451_1;
7634                                      }                                      }
7635                                  }                                  }
7636                              } else { /* constant data */                              } else { // constant data
7637                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
7638                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
7639                                          const register double B_0 = B_p[INDEX3(k,0,m, numEq, 3)];                                          const register double B_0 = B_p[INDEX3(k,0,m, numEq, 3)];
# Line 8357  void Brick::assemblePDESystem(Paso_Syste Line 8322  void Brick::assemblePDESystem(Paso_Syste
8322                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp388_1 + tmp389_1 + tmp390_1 + tmp391_1 + tmp396_1 + tmp397_1 + tmp448_1 + tmp449_1 + tmp450_1 + tmp451_1;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp388_1 + tmp389_1 + tmp390_1 + tmp391_1 + tmp396_1 + tmp397_1 + tmp448_1 + tmp449_1 + tmp450_1 + tmp451_1;
8323                                      }                                      }
8324                                  }                                  }
8325                              } else { /* constant data */                              } else { // constant data
8326                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8327                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8328                                          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 8640  void Brick::assemblePDESystem(Paso_Syste Line 8605  void Brick::assemblePDESystem(Paso_Syste
8605                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp71_1 + tmp72_1;                                          EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp71_1 + tmp72_1;
8606                                      }                                      }
8607                                   }                                   }
8608                              } else { /* constant data */                              } else { // constant data
8609                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8610                                      for (index_t m=0; m<numComp; m++) {                                      for (index_t m=0; m<numComp; m++) {
8611                                          const register double D_0 = D_p[INDEX2(k, m, numEq)];                                          const register double D_0 = D_p[INDEX2(k, m, numEq)];
# Line 8835  void Brick::assemblePDESystem(Paso_Syste Line 8800  void Brick::assemblePDESystem(Paso_Syste
8800                                      EM_F[INDEX2(k,6,numEq)]+=tmp28_1 + tmp2_1 + tmp43_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1 + tmp52_1 + tmp53_1;                                      EM_F[INDEX2(k,6,numEq)]+=tmp28_1 + tmp2_1 + tmp43_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1 + tmp52_1 + tmp53_1;
8801                                      EM_F[INDEX2(k,7,numEq)]+=tmp12_1 + tmp21_1 + tmp37_1 + tmp54_1 + tmp55_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;                                      EM_F[INDEX2(k,7,numEq)]+=tmp12_1 + tmp21_1 + tmp37_1 + tmp54_1 + tmp55_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;
8802                                  }                                  }
8803                              } else { /* constant data */                              } else { // constant data
8804                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8805                                      const register double X_0 = X_p[INDEX2(k, 0, numEq)];                                      const register double X_0 = X_p[INDEX2(k, 0, numEq)];
8806                                      const register double X_1 = X_p[INDEX2(k, 1, numEq)];                                      const register double X_1 = X_p[INDEX2(k, 1, numEq)];
# Line 8922  void Brick::assemblePDESystem(Paso_Syste Line 8887  void Brick::assemblePDESystem(Paso_Syste
8887                                      EM_F[INDEX2(k,6,numEq)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1;                                      EM_F[INDEX2(k,6,numEq)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1;
8888                                      EM_F[INDEX2(k,7,numEq)]+=tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;                                      EM_F[INDEX2(k,7,numEq)]+=tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;
8889                                  }                                  }
8890                              } else { /* constant data */                              } else { // constant data
8891                                  for (index_t k=0; k<numEq; k++) {                                  for (index_t k=0; k<numEq; k++) {
8892                                      const register double Y_0 = Y_p[k];                                      const register double Y_0 = Y_p[k];
8893                                      const register double tmp0_1 = Y_0*w183;                                      const register double tmp0_1 = Y_0*w183;
# Line 8937  void Brick::assemblePDESystem(Paso_Syste Line 8902  void Brick::assemblePDESystem(Paso_Syste
8902                                  }                                  }
8903                              }                              }
8904                          }                          }
                         /* GENERATOR SNIP_PDE_SYSTEM BOTTOM */  
8905    
8906                          // 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)
8907                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
# Line 8951  void Brick::assemblePDESystem(Paso_Syste Line 8915  void Brick::assemblePDESystem(Paso_Syste
8915                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);
8916                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);
8917                          if (add_EM_F) {                          if (add_EM_F) {
                             //cout << "-- AddtoRHS -- " << endl;  
8918                              double *F_p=rhs.getSampleDataRW(0);                              double *F_p=rhs.getSampleDataRW(0);
8919                              for (index_t i=0; i<rowIndex.size(); i++) {                              for (index_t i=0; i<rowIndex.size(); i++) {
8920                                  if (rowIndex[i]<getNumDOF()) {                                  if (rowIndex[i]<getNumDOF()) {
8921                                      for (index_t eq=0; eq<numEq; eq++) {                                      for (index_t eq=0; eq<numEq; eq++) {
8922                                          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;  
8923                                      }                                      }
8924                                  }                                  }
8925                              }                              }
                             //cout << "---"<<endl;  
8926                          }                          }
8927                          if (add_EM_S) {                          if (add_EM_S) {
                             //cout << "-- AddtoSystem -- " << endl;  
8928                              addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);                              addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);
8929                          }                          }
8930                      } // end k0 loop                      } // end k0 loop
# Line 8992  void Brick::assemblePDESystemReduced(Pas Line 8952  void Brick::assemblePDESystemReduced(Pas
8952          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
8953      }      }
8954    
8955      /* GENERATOR SNIP_PDE_SYSTEM_REDUCED_PRE TOP */      const double w0 = .0625*h1*h2/h0;
8956      const double w0 = 0.0625*h1*h2/h0;      const double w1 = .0625*h2;
8957      const double w1 = 0.0625*h2;      const double w2 = -.0625*h1;
8958      const double w10 = -0.0625*h0*h2/h1;      const double w3 = .0625*h0*h2/h1;
8959      const double w11 = 0.0625*h0*h1/h2;      const double w4 = -.0625*h0;
8960      const double w12 = 0.03125*h1*h2;      const double w5 = .0625*h1;
8961      const double w13 = 0.03125*h0*h2;      const double w6 = .0625*h0;
8962      const double w14 = 0.03125*h0*h1;      const double w7 = -.0625*h0*h1/h2;
8963      const double w15 = -0.03125*h1*h2;      const double w8 = -.0625*h1*h2/h0;
8964      const double w16 = -0.03125*h0*h2;      const double w9 = -.0625*h2;
8965      const double w17 = -0.03125*h0*h1;      const double w10 = -.0625*h0*h2/h1;
8966      const double w18 = 0.015625*h0*h1*h2;      const double w11 = .0625*h0*h1/h2;
8967      const double w19 = -0.25*h1*h2;      const double w12 = .03125*h1*h2;
8968      const double w2 = -0.0625*h1;      const double w13 = .03125*h0*h2;
8969      const double w20 = -0.25*h0*h2;      const double w14 = .03125*h0*h1;
8970      const double w21 = -0.25*h0*h1;      const double w15 = -.03125*h1*h2;
8971      const double w22 = 0.25*h1*h2;      const double w16 = -.03125*h0*h2;
8972      const double w23 = 0.25*h0*h2;      const double w17 = -.03125*h0*h1;
8973      const double w24 = 0.25*h0*h1;      const double w18 = .015625*h0*h1*h2;
8974      const double w25 = 0.125*h0*h1*h2;      const double w19 = -.25*h1*h2;
8975      const double w3 = 0.0625*h0*h2/h1;      const double w20 = -.25*h0*h2;
8976      const double w4 = -0.0625*h0;      const double w21 = -.25*h0*h1;
8977      const double w5 = 0.0625*h1;      const double w22 = .25*h1*h2;
8978      const double w6 = 0.0625*h0;      const double w23 = .25*h0*h2;
8979      const double w7 = -0.0625*h0*h1/h2;      const double w24 = .25*h0*h1;
8980      const double w8 = -0.0625*h1*h2/h0;      const double w25 = .125*h0*h1*h2;
     const double w9 = -0.0625*h2;  
     /* GENERATOR SNIP_PDE_SYSTEM_REDUCED_PRE BOTTOM */  
8981    
8982      rhs.requireWrite();      rhs.requireWrite();
8983  #pragma omp parallel  #pragma omp parallel
# Line 9034  void Brick::assemblePDESystemReduced(Pas Line 8992  void Brick::assemblePDESystemReduced(Pas
8992                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
8993                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
8994                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;
                         /* GENERATOR SNIP_PDE_SYSTEM_REDUCED TOP */  
8995                          ///////////////                          ///////////////
8996                          // process A //                          // process A //
8997                          ///////////////                          ///////////////
# Line 9436  void Brick::assemblePDESystemReduced(Pas Line 9393  void Brick::assemblePDESystemReduced(Pas
9393                                  EM_F[INDEX2(k,7,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,7,numEq)]+=tmp0_1;
9394                              }                              }
9395                          }                          }
                         /* GENERATOR SNIP_PDE_SYSTEM_REDUCED BOTTOM */  
9396    
9397                          // 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)
9398                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
# Line 9450  void Brick::assemblePDESystemReduced(Pas Line 9406  void Brick::assemblePDESystemReduced(Pas
9406                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);
9407                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);                          rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);
9408                          if (add_EM_F) {                          if (add_EM_F) {
                             //cout << "-- AddtoRHS -- " << endl;  
9409                              double *F_p=rhs.getSampleDataRW(0);                              double *F_p=rhs.getSampleDataRW(0);
9410                              for (index_t i=0; i<rowIndex.size(); i++) {                              for (index_t i=0; i<rowIndex.size(); i++) {
9411                                  if (rowIndex[i]<getNumDOF()) {                                  if (rowIndex[i]<getNumDOF()) {
9412                                      for (index_t eq=0; eq<numEq; eq++) {                                      for (index_t eq=0; eq<numEq; eq++) {
9413                                          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;  
9414                                      }                                      }
9415                                  }                                  }
9416                              }                              }
                             //cout << "---"<<endl;  
9417                          }                          }
9418                          if (add_EM_S) {                          if (add_EM_S) {
                             //cout << "-- AddtoSystem -- " << endl;  
9419                              addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);                              addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);
9420                          }                          }
9421                      } // end k0 loop                      } // end k0 loop

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

  ViewVC Help
Powered by ViewVC 1.1.26