Actual source code: FEMProblem.hh
petsc-3.4.2 2013-07-02

      1: #ifndef __FEMProblem

      4: /*

      6:   Framework for solving a FEM problem using sieve.

      8:   The Discretization objects are WAY too embedded into the way things are done; we will have to create ways of

     10: This includes derived types doing what indicesExcluded does for all things marked with a boundary marker.

     12: */

     14: #include <sieve/Mesh.hh>

     16: namespace ALE {
     17:   namespace Problem {

     19:     /*
     20:       This class defines a basic interface to a subproblem; all data the type needs will be set at initialization and be
     21:       members of the derived types

     23:       The recommended use for a given problem is to define a subproblem class for doing the data initialization for the
     24:       form, as well as the postprocessing one for doing the boundary handling and setting, and one for dealing with the
     25:       process of the solve.  Work can be broken up as needed; however for the sake of simplicity one should probably
     26:       have various forms in various subproblems to assemble.  This could of course include facet integrals or problems
     27:       over part of the domain.  Look at the examples provided in UFCProblem, which use the UFC form compiler to assemble
     28:       subproblems of a form.

     30:      */

     32:     class SubProblem : public ParallelObject {
     33:     public:

     35:       typedef std::string name_type;

     37:       SubProblem(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {};
     38:       virtual ~SubProblem() {};

     40:     };

     42:     /*
     43:       Below is a helper derived class of subproblem; the persistently limiting and annoying discretization object in
     44:       Sieve must be excised; Here are subproblem objects containing the helper functions that are presently in the
     45:       monster called PETSC_MESH_TYPE, but with references to the discretization object stolen.  virtual functions to be
     46:       filled in in further derived classes are marked.  This is meant to be a stepping stone towards generalized use of
     47:       problem creation objects with multiple fields and interesting boundary conditions.  This type can be branched into

     49:     */

     51:     //creates a discretization-like thing for this particular case; from this we can use the helper functions in the
     52:     //generalformsubproblem to set up the discretization across the mesh as one might expect in the case of a

     55: #if 0

     57:     //we really can't use this for finite elements with UFC

     59:     class GeneralCell : public ParallelObject {
     60:     private:
     61:       int _embedded_dimension;           //the fiberdimension of the coordinate section
     62:       int closureSize;                   //the size of the closure
     63:       std::map<int, int> _closure_order; //all that really matters; assume interpolated as the order should be the same
     64:       int _num_vertices;                 //necessary to allocate the coordinate array.
     65:       double * _coordinates;             //the coordinate array in order based upon the local topology

     67:       GeneralCell() {
     68:         _num_vertices = 0;
     69:         _coordinates = NULL;
     70:       }

     72:       GeneralCell(int embedded_dimension, int num_vertices) {
     73:         _embedded_dimension = embedded_dimension;
     74:         _num_vertices = _num_vertices;
     75:         _coordinates = new double[_embedded_dimension*_num_vertices];
     76:       }

     78:       virtual ~GeneralCell() {
     79:         if (_coordinates)
     80:           delete _coordinates;
     81:       }

     83:       virtual void setMeshCell(Obj<PETSC_MESH_TYPE> mesh, PETSC_MESH_TYPE::point_type cell) {
     84:         throw Exception("GeneralCell->setCell(): Unimplemented base class");
     85:         return;
     86:       }
     87:       virtual void setClosureOrder(int subcell, int map) {
     88:         _closure_order[subcell] = map;
     89:       }
     90:       virtual int getClosureOrder(int subcell) {
     91:         return _closure_order[subcell];
     92:       }
     93:     };

     95: #endif

     97:     class GeneralBoundaryCondition : ParallelObject {
     98:     protected:

    100:       std::string     _labelName;
    101:       int             _marker;

    103:     public:
    104:       GeneralBoundaryCondition(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {};
    105:       virtual ~GeneralBoundaryCondition() {};
    106:     public:
    107:       virtual const std::string& getLabelName() const {return this->_labelName;};
    108:       virtual void setLabelName(const std::string& name) {this->_labelName = name;};
    109:       virtual int getMarker() const {return this->_marker;};
    110:       virtual void setMarker(const int marker) {this->_marker = marker;};

    112:     public:

    114:       virtual double integrateDual(unsigned int dof) {
    115:         //when implementing this one should have some notion of what the dof number is built into the BC; this would constitute having some cell or something known
    116:         //by the object such that that cell can be set and integrated within.
    117:         throw Exception("GeneralBoundaryCondition->integrateDual: Nonimplemented base-class version called.");
    118:         return 3.;
    119:       };

    121:       virtual void setReorder(int * reorder) {
    122:         throw Exception("GeneralBoundaryCondition->setReorder(): Unimplemented base class version called.");
    123:       }

    125:       virtual const int * getReorder() {
    126:         throw Exception("GeneralBoundaryCondition->getReorder(): Unimplemented base class version called.");
    127:         return NULL;
    128:       }

    130:     };

    132:     /*
    133:       Include at least counts of all the part of the triple, as well as all the information for the cell.
    134:      */

    136: #if 0

    138:     class GeneralFiniteElement : public ParallelObject {
    139:     private:
    140:     public:
    141:       virtual double integrateDual(unsigned int dof) {
    142:         //evaluate a degree of freedom
    143:         return 0.;
    144:       }
    145:       virtual int closureIndex(unsigned int dof) {
    146:         return 0;
    147:       }
    148:       virtual int dataIndex(unsigned int dof) {
    149:         return 0;
    150:       }
    151:     };

    153: #endif

    155:     //we almost, ALMOST need an overall view of a local reference topology for this kind of stuff (and the boundary conditions).

    157:     class GeneralIntegral : public ParallelObject {
    158:     public:
    159:       typedef std::string name_type;
    160:     private:
    161:       //the integral should only apply to the labeled points in some way, probably over the closure(support(point));

    164:       //the integrals will have some set of cells or cell-like objects they act on; this includes potentially internal faces.
    165:       //But, for cell integrals these should be set to "height" and "0"

    167:       name_type _label_name;
    168:       int _label_marker;     //We could use this if there's only a certain label and marker that the given integral applies to

    170:       int _num_coefficients; //for the linear forms

    172:       int _space_dimension;       //the number of DoFs we're dealing with here.
    173:       int _tensor_rank;           //the tensor rank; it BETTER be one or two.
    174:       int _topological_dimension; //the topological dimension of the given mesh item -- tells us if it's a cell or facet integral

    176:       int * _closure2data;  //if there is some API-level data storage, this maps the unknowns for the WHOLE CLOSURE onto the unknowns for the API

    178:     public:

    180:       GeneralIntegral (MPI_Comm comm, int debug = 0) : ParallelObject(comm, debug) {
    181:         _tensor_rank = 0;
    182:         _space_dimension = 0;
    183:         _topological_dimension = 0;
    184:         _label_marker = 0;
    185:       }

    187:       GeneralIntegral (MPI_Comm comm, name_type labelname, int labelmarker, int nCoefficients = 0, int debug = 0) : ParallelObject(comm, debug) {
    188:         _num_coefficients = nCoefficients;
    189:         _label_name = labelname;
    190:         _label_marker = labelmarker;
    191:       }

    193:       virtual ~GeneralIntegral () {};

    196:       int getNumCoefficients() {
    197:         return _num_coefficients;
    198:       }

    200:       void setNumCoefficients(int nc) {
    201:         _num_coefficients = nc;
    202:       }

    205:       name_type getLabelName() {
    206:         return _label_name;
    207:       }

    209:       int getLabelMarker() {
    210:         return _label_marker;
    211:       }

    213:       void setLabelName(name_type newName) {
    214:         _label_name = newName;
    215:       }

    217:       void setLabelMarker(int newMarker) {
    218:         _label_marker = newMarker;
    219:       }

    222:       virtual void setSpaceDimension(int dim) {
    223:         _space_dimension = dim;
    224:       }

    226:       virtual int getSpaceDimension() {
    227:         return _space_dimension;
    228:       }

    230:       virtual void setTensorRank(int rank) {
    231:         _tensor_rank = rank;
    232:       }

    234:       virtual int getTensorRank() {
    235:         return _tensor_rank;
    236:       }

    240:       virtual const int * getReorder() {
    241:         return _closure2data;
    242:       }

    244:       virtual void setReorder(int * reorder) {
    245:         _closure2data = reorder;
    246:       }

    248:       //use the UFC lingo here
    249:       //have some notion of the cell initialized and in-state in the eventual implementation.
    250:       virtual void tabulateTensor(double * tensor, const double * coefficients = NULL) {
    251:         throw Exception("GeneralIntegral->tabulateTensor: Nonimplemented Base class version called.");
    252:         return;
    253:       }
    254:     };

    256:     class GeneralDiscretization : ParallelObject { //this should largely resemble the old form of the discretizations, only with specifics of evaluation to FIAT removed.
    257:       typedef std::map<std::string, Obj<GeneralBoundaryCondition> > boundaryConditions_type;
    258:     protected:
    259:       boundaryConditions_type _boundaryConditions;
    260:       Obj<GeneralBoundaryCondition> _exactSolution;
    261:       std::map<int,int> _dim2dof; //not good enough for tensor product assembly.... however we can generalize this into some sort of "getClosureItemDofs" or something per cell
    262:       std::map<int,int> _dim2class;
    263:       int           _quadSize;
    264:       int           _basisSize;
    265:       const int *         _indices;
    266:       std::map<int, const int *> _exclusionIndices;
    267:       const int * _closure2data; //local index reordering

    269:     public:

    271:       typedef std::set<std::string> names_type;

    273:       GeneralDiscretization(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _quadSize(0), _basisSize(0), _indices(NULL) {};
    274:       virtual ~GeneralDiscretization() {
    275:         if (this->_indices) {delete [] this->_indices;}
    276:         for(std::map<int, const int *>::iterator i_iter = _exclusionIndices.begin(); i_iter != _exclusionIndices.end(); ++i_iter) {
    277:           delete [] i_iter->second;
    278:         }
    279:         if (_closure2data) {
    280:           delete _closure2data;
    281:         }
    282:       };
    283:     public:
    284:       virtual const bool hasBoundaryCondition() {return (this->_boundaryConditions.find("default") != this->_boundaryConditions.end());};
    285:       virtual Obj<GeneralBoundaryCondition> getBoundaryCondition() {return this->getBoundaryCondition("default");};
    286:       virtual void setBoundaryCondition(const Obj<GeneralBoundaryCondition>& boundaryCondition) {this->setBoundaryCondition("default", boundaryCondition);};
    287:       virtual Obj<GeneralBoundaryCondition> getBoundaryCondition(const std::string& name) {return this->_boundaryConditions[name];};
    288:       virtual void setBoundaryCondition(const std::string& name, const Obj<GeneralBoundaryCondition>& boundaryCondition) {this->_boundaryConditions[name] = boundaryCondition;};
    289:       virtual names_type getBoundaryConditions() const {
    290:         Obj<names_type> names = names_type();
    291:         for(boundaryConditions_type::const_iterator d_iter = this->_boundaryConditions.begin(); d_iter != this->_boundaryConditions.end(); ++d_iter) {
    292:           names->insert(d_iter->first);
    293:         }
    294:         return names;
    295:       };

    297:       //eh?
    298:       //virtual const Obj<BoundaryCondition>& getExactSolution() {return this->_exactSolution;};
    299:       //virtual void setExactSolution(const Obj<GeneralBoundaryCondition>& exactSolution) {this->_exactSolution = exactSolution;};
    300:       virtual const int     getQuadratureSize() {return this->_quadSize;};
    301:       virtual void          setQuadratureSize(const int size) {this->_quadSize = size;};
    302:       virtual const int     getBasisSize() {return this->_basisSize;};
    303:       virtual void          setBasisSize(const int size) {this->_basisSize = size;};

    305:       //eh, until they get this together in UFC keep these around.

    307:       virtual int           getNumDof(const int dim) {return this->_dim2dof[dim];};
    308:       virtual void          setNumDof(const int dim, const int numDof) {this->_dim2dof[dim] = numDof;};
    309:       virtual int           getDofClass(const int dim) {return this->_dim2class[dim];};
    310:       virtual void          setDofClass(const int dim, const int dofClass) {this->_dim2class[dim] = dofClass;};
    311:     public:

    313:       /*
    314:
    315:       Functions for interacting with external libraries for handling finite element assembly that might have different cell layout.

    317:        */

    319:       virtual void createReorder() {
    320:         throw Exception("GeneralDiscretization->createReorderings: Unimplemented base function");
    321:         return;
    322:       }

    324:       virtual const int * getReorder() {
    325:         return _closure2data;
    326:       }

    328:       virtual void setReorder(int * reorder) {
    329:         _closure2data = reorder;
    330:       }

    332:       //Yeah... not messing with this part.

    334:       virtual const int *getIndices() {return this->_indices;};
    335:       virtual const int *getIndices(const int marker) {

    337:         if (!marker) return this->getIndices();
    338:         return this->_exclusionIndices[marker];
    339:       };
    340:       virtual void       setIndices(const int *indices) {this->_indices = indices;};
    341:       virtual void       setIndices(const int *indices, const int marker) {
    342:         if (!marker) this->_indices = indices;
    343:         this->_exclusionIndices[marker] = indices;
    344:       };

    346:       //return the size of the space
    347:       virtual int size() {
    348:         throw Exception("GeneralDiscretization->size(): Nonimplemented base class function called.");
    349:         return 0;
    350:       }
    351:       virtual double evaluateRHS(int dof) {
    352:         throw Exception("GeneralDiscretization->evaluateRHS: Nonimplemented base class function called.");
    353:       }
    354:     };

    356:     //The GenericFormSubProblem shown here should basically contain the whole problem as it does the discretization-like setup, which might break
    357:     //if you define multiple ones right now.  Set the multiple different spaces using different discretizations.

    359:     class GenericFormSubProblem : SubProblem { //takes information from some abstract notion of a problem and sets up the thing.
    360:     public:
    361:       typedef std::map<std::string, Obj<GeneralDiscretization> > discretizations_type;
    362:       typedef std::map<std::string, Obj<GeneralIntegral> > integral_type;
    363:       typedef std::set<std::string> names_type;

    365:       GenericFormSubProblem(MPI_Comm comm, int debug = 1) : SubProblem(comm, debug) {};

    367:       ~GenericFormSubProblem(){};

    369:     private:

    371:       name_type _solutionSectionName;
    372:       name_type _forcingSectionName;
    373:       discretizations_type _discretizations;
    374:       integral_type _integrals;
    375:       Obj<GeneralBoundaryCondition> _exactSolution; //evaluates a function over all unknowns on the form containing an exact solution.  Per discretization later.
    376:       int * _closure2data;                          //a mapping from closure indices to data indices for the overall element

    379:       //helper functions, stolen directly from Mesh.hh with slight modifications to remove FIAT dependence and instead use stuff from GeneralDiscretization

    381:     public:

    383:       const Obj<GeneralDiscretization>& getDiscretization() {return this->getDiscretization("default");};
    384:       const Obj<GeneralDiscretization>& getDiscretization(const std::string& name) {return this->_discretizations[name];};

    386:       void setDiscretization(const Obj<GeneralDiscretization>& disc) {this->setDiscretization("default", disc);};
    387:       void setDiscretization(const std::string& name, const Obj<GeneralDiscretization>& disc) {this->_discretizations[name] = disc;};

    389:       const Obj<GeneralIntegral>& getIntegral() {return this->getIntegral("default");};
    390:       const Obj<GeneralIntegral>& getIntegral(const std::string& name) {return this->_integrals[name];};

    392:       void setIntegral(const Obj<GeneralIntegral>& integral) {this->setIntegral("default", integral);};
    393:       void setIntegral(const std::string& name, const Obj<GeneralIntegral> integral) {this->_integrals[name] = integral;};

    395:       //FAIL! this won't be used.
    396:       void setExactSolution(Obj<GeneralBoundaryCondition> exactsol) {
    397:         _exactSolution = exactsol;
    398:       }

    400:       Obj<GeneralBoundaryCondition> getExactSolution(){return _exactSolution;};

    402:       Obj<names_type> getDiscretizations() const {
    403:         Obj<names_type> names = names_type();
    404:
    405:         for(discretizations_type::const_iterator d_iter = this->_discretizations.begin(); d_iter != this->_discretizations.end(); ++d_iter) {
    406:           names->insert(d_iter->first);
    407:         }
    408:         return names;
    409:       };

    411:       Obj<names_type> getIntegrals() const {
    412:         Obj<names_type> names = names_type();
    413:         for (integral_type::const_iterator i_iter = this->_integrals.begin(); i_iter != this->_integrals.end(); ++i_iter) {
    414:           names->insert(i_iter->first);
    415:         }
    416:         return names;
    417:       }

    419:       //TODO: Implement a cell interface class towards flexible handling of reordering between various parts of this thing.

    421:       virtual void setCell(Obj<PETSC_MESH_TYPE> mesh, PETSC_MESH_TYPE::point_type c) const { //implementations should set some state of the derived type to the present cell
    422:         throw Exception("Using the unimplemented base class version of setCell in GeneralDiscretization.");
    423:         return;
    424:       };

    426:       virtual int localSpaceDimension() {
    427:         throw Exception("GeneralDiscretization->localSpaceDimension(): using the unimplemented base class version.");
    428:         return 0;
    429:       }

    431:       /*
    432:         Functions handling the creation and application of reordering from element libraries and sieve.
    433:        */

    435:       virtual void createReorder() {
    436:         throw Exception("GeneralFormSubProblem->buildOrderings(): nonimplemented base function");
    437:         return;
    438:       }

    440:       virtual const int * getReorder() {
    441:         return _closure2data;
    442:       }

    444:       virtual void setReorder(int * order) {
    445:         _closure2data = order;
    446:       }

    448:       /*
    449:         Functions handling the mesh data layout from the overall subproblem
    450:        */

    452:       int setFiberDimensions(const Obj<PETSC_MESH_TYPE> mesh, const Obj<PETSC_MESH_TYPE::real_section_type> s, const Obj<names_type>& discs, names_type& bcLabels) {
    453:         const int debug  = this->debug();
    454:         int       maxDof = 0;
    455:
    456:         for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
    457:           s->addSpace();
    458:         }
    459:         for(int d = 0; d <= mesh->getDimension(); ++d) {
    460:           int numDof = 0;
    461:           int f      = 0;
    462:
    463:           for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
    464:             const Obj<GeneralDiscretization>& disc = this->getDiscretization(*f_iter);
    465:             const int                       sDof = disc->getNumDof(d);
    466:
    467:             numDof += sDof;
    468:             if (sDof) s->setFiberDimension(mesh->depthStratum(d), sDof, f);
    469:           }
    470:           if (numDof) s->setFiberDimension(mesh->depthStratum(d), numDof);
    471:           maxDof = std::max(maxDof, numDof);
    472:         }
    473:         // Process exclusions
    474:         typedef ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> Visitor;
    475:         int f = 0;
    476:
    477:         for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
    478:           const Obj<GeneralDiscretization>& disc      = this->getDiscretization(*f_iter);
    479:           std::string                     labelName = "exclude-"+*f_iter;
    480:           std::set<PETSC_MESH_TYPE::point_type>            seen;
    481:           Visitor pV((int) pow(mesh->getSieve()->getMaxConeSize(), mesh->depth()), true);
    482:
    483:           if (mesh->hasLabel(labelName)) {
    484:             const Obj<PETSC_MESH_TYPE::label_type>&         label     = mesh->getLabel(labelName);
    485:             const Obj<PETSC_MESH_TYPE::label_sequence>&     exclusion = mesh->getLabelStratum(labelName, 1);
    486:             const PETSC_MESH_TYPE::label_sequence::iterator end       = exclusion->end();
    487:             if (debug > 1) {label->view(labelName.c_str());}
    488:
    489:             for(PETSC_MESH_TYPE::label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
    490:               ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mesh->getSieve(), *e_iter, pV);
    491:               const Visitor::point_type *oPoints = pV.getPoints();
    492:               const int                  oSize   = pV.getSize();
    493:
    494:               for(int cl = 0; cl < oSize; ++cl) {
    495:                 if (seen.find(oPoints[cl]) != seen.end()) continue;
    496:                 if (mesh->getValue(label, oPoints[cl]) == 1) {
    497:                   seen.insert(oPoints[cl]);
    498:                   s->setFiberDimension(oPoints[cl], 0, f);
    499:                   s->addFiberDimension(oPoints[cl], -disc->getNumDof(mesh->depth(oPoints[cl])));
    500:                   if (debug > 1) {std::cout << "  point: " << oPoints[cl] << " dim: " << disc->getNumDof(mesh->depth(oPoints[cl])) << std::endl;}
    501:                 }
    502:               }
    503:               pV.clear();
    504:             }
    505:           }
    506:         }
    507:         // Process constraints
    508:         f = 0;
    509:         for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
    510:           const Obj<GeneralDiscretization>&    disc        = this->getDiscretization(*f_iter);
    511:           const Obj<std::set<std::string> >  bcs         = disc->getBoundaryConditions();
    512:           std::string                        excludeName = "exclude-"+*f_iter;
    513:
    514:           for(std::set<std::string>::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter) {
    515:             const Obj<GeneralBoundaryCondition>& bc       = disc->getBoundaryCondition(*bc_iter);
    516:             const Obj<PETSC_MESH_TYPE::label_sequence>&         boundary = mesh->getLabelStratum(bc->getLabelName(), bc->getMarker());
    517:
    518:             bcLabels.insert(bc->getLabelName());
    519:             if (mesh->hasLabel(excludeName)) {
    520:               const Obj<PETSC_MESH_TYPE::label_type>& label = mesh->getLabel(excludeName);
    521:
    522:               for(PETSC_MESH_TYPE::label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
    523:                 if (!mesh->getValue(label, *e_iter)) {
    524:                   const int numDof = disc->getNumDof(mesh->depth(*e_iter));
    525:
    526:                   if (numDof) s->addConstraintDimension(*e_iter, numDof);
    527:                   if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
    528:                 }
    529:               }
    530:             } else {
    531:               for(PETSC_MESH_TYPE::label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
    532:                 const int numDof = disc->getNumDof(mesh->depth(*e_iter));
    533:
    534:                 if (numDof) s->addConstraintDimension(*e_iter, numDof);
    535:                 if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
    536:               }
    537:             }
    538:           }
    539:         }
    540:         return maxDof;
    541:       };
    542:       void calculateIndices(Obj<PETSC_MESH_TYPE> mesh) {
    543:         typedef ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> Visitor;
    544:         // Should have an iterator over the whole tree
    545:         Obj<names_type> discs = this->getDiscretizations();
    546:         const int       debug = this->debug();
    547:         std::map<std::string, std::pair<int, int*> > indices;
    548:
    549:         for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
    550:           const Obj<GeneralDiscretization>& disc = this->getDiscretization(*d_iter);
    551:           indices[*d_iter] = std::pair<int, int*>(0, new int[disc->size()]);  //size isn't a function of the mesh it's a function of the local function space
    552:           disc->setIndices(indices[*d_iter].second);
    553:         }
    554:         const Obj<PETSC_MESH_TYPE::label_sequence>& cells   = mesh->heightStratum(0);
    555:         Visitor pV((int) pow(mesh->getSieve()->getMaxConeSize(), mesh->depth())+1, true);
    556:         ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mesh->getSieve(), *cells->begin(), pV);
    557:         const Visitor::point_type *oPoints = pV.getPoints();
    558:         const int                  oSize   = pV.getSize();
    559:         int                        offset  = 0;
    560:
    561:         if (debug > 1) {std::cout << "Closure for first element" << std::endl;}
    562:         for(int cl = 0; cl < oSize; ++cl) {
    563:           const int dim = mesh->depth(oPoints[cl]);
    564:
    565:           if (debug > 1) {std::cout << "  point " << oPoints[cl] << " depth " << dim << std::endl;}
    566:           for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
    567:             const Obj<GeneralDiscretization>& disc = this->getDiscretization(*d_iter);
    568:             const int                  num  = disc->getNumDof(dim);
    569:
    570:             //if (debug > 1) {std::cout << "    disc " << disc->getName() << " numDof " << num << std::endl;}
    571:             for(int o = 0; o < num; ++o) {
    572:               indices[*d_iter].second[indices[*d_iter].first++] = offset++;
    573:             }
    574:           }
    575:         }
    576:         pV.clear();
    577:         if (debug > 1) {
    578:           for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
    579:             //const Obj<GeneralDiscretization>& disc = this->getDiscretization(*d_iter);
    580:
    581:             //std::cout << "Discretization " << disc->getName() << " indices:";
    582:             for(int i = 0; i < indices[*d_iter].first; ++i) {
    583:               std::cout << " " << indices[*d_iter].second[i];
    584:             }
    585:             std::cout << std::endl;
    586:           }
    587:         }
    588:       };

    590:       void calculateIndicesExcluded(const Obj<PETSC_MESH_TYPE> mesh, const Obj<PETSC_MESH_TYPE::real_section_type>& s, const Obj<names_type>& discs) {
    591:         typedef ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> Visitor;
    592:         typedef std::map<std::string, std::pair<int, indexSet> > indices_type;
    593:         const Obj<PETSC_MESH_TYPE::label_type>& indexLabel = mesh->createLabel("cellExclusion");
    594:         const int debug  = this->debug();
    595:         int       marker = 0;
    596:         std::map<indices_type, int> indexMap;
    597:         indices_type                indices;
    598:         Visitor pV((int) pow(mesh->getSieve()->getMaxConeSize(), mesh->depth())+1, true);
    599:
    600:         for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
    601:           const Obj<GeneralDiscretization>& disc = this->getDiscretization(*d_iter);
    602:           const int                  size = disc->size();
    603:
    604:           indices[*d_iter].second.resize(size);
    605:         }
    606:         const names_type::const_iterator dBegin = discs->begin();
    607:         const names_type::const_iterator dEnd   = discs->end();
    608:         std::set<PETSC_MESH_TYPE::point_type> seen;
    609:         int f = 0;
    610:
    611:         for(names_type::const_iterator f_iter = dBegin; f_iter != dEnd; ++f_iter, ++f) {
    612:           std::string labelName = "exclude-"+*f_iter;
    613:
    614:           if (mesh->hasLabel(labelName)) {
    615:             const Obj<PETSC_MESH_TYPE::label_sequence>&     exclusion = mesh->getLabelStratum(labelName, 1);
    616:             const PETSC_MESH_TYPE::label_sequence::iterator end       = exclusion->end();
    617:
    618:             if (debug > 1) {std::cout << "Processing exclusion " << labelName << std::endl;}
    619:             for(PETSC_MESH_TYPE::label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
    620:               if (mesh->height(*e_iter)) continue;
    621:               ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mesh->getSieve(), *e_iter, pV);
    622:               const Visitor::point_type *oPoints = pV.getPoints();
    623:               const int                  oSize   = pV.getSize();
    624:               int                        offset  = 0;
    625:
    626:               if (debug > 1) {std::cout << "  Closure for cell " << *e_iter << std::endl;}
    627:               for(int cl = 0; cl < oSize; ++cl) {
    628:                 int g = 0;
    629:
    630:                 if (debug > 1) {std::cout << "    point " << oPoints[cl] << std::endl;}
    631:                 for(names_type::const_iterator g_iter = dBegin; g_iter != dEnd; ++g_iter, ++g) {
    632:                   const int fDim = s->getFiberDimension(oPoints[cl], g);
    633:
    634:                   if (debug > 1) {std::cout << "      disc " << *g_iter << " numDof " << fDim << std::endl;}
    635:                   for(int d = 0; d < fDim; ++d) {
    636:                     indices[*g_iter].second[indices[*g_iter].first++] = offset++;
    637:                   }
    638:                 }
    639:               }
    640:               pV.clear();
    641:               const std::map<indices_type, int>::iterator entry = indexMap.find(indices);
    642:
    643:               if (debug > 1) {
    644:                 for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
    645:                   for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
    646:                     std::cout << "Discretization (" << i_iter->second << ") " << *g_iter << " indices:";
    647:                     for(int i = 0; i < ((indices_type) i_iter->first)[*g_iter].first; ++i) {
    648:                       std::cout << " " << ((indices_type) i_iter->first)[*g_iter].second[i];
    649:                     }
    650:                     std::cout << std::endl;
    651:                   }
    652:                   std::cout << "Comparison: " << (indices == i_iter->first) << std::endl;
    653:                 }
    654:                 for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
    655:                   std::cout << "Discretization " << *g_iter << " indices:";
    656:                   for(int i = 0; i < indices[*g_iter].first; ++i) {
    657:                     std::cout << " " << indices[*g_iter].second[i];
    658:                   }
    659:                   std::cout << std::endl;
    660:                 }
    661:               }
    662:               if (entry != indexMap.end()) {
    663:                 mesh->setValue(indexLabel, *e_iter, entry->second);
    664:                 if (debug > 1) {std::cout << "  Found existing indices with marker " << entry->second << std::endl;}
    665:               } else {
    666:                 indexMap[indices] = ++marker;
    667:                 mesh->setValue(indexLabel, *e_iter, marker);
    668:                 if (debug > 1) {std::cout << "  Created new indices with marker " << marker << std::endl;}
    669:               }
    670:               for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
    671:                 indices[*g_iter].first  = 0;
    672:                 for(unsigned int i = 0; i < indices[*g_iter].second.size(); ++i) indices[*g_iter].second[i] = 0;
    673:               }
    674:             }
    675:           }
    676:         }
    677:         if (debug > 1) {indexLabel->view("cellExclusion");}
    678:         for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
    679:           if (debug > 1) {std::cout << "Setting indices for marker " << i_iter->second << std::endl;}
    680:           for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
    681:             const Obj<GeneralDiscretization>& disc = this->getDiscretization(*g_iter);
    682:             const indexSet  indSet   = ((indices_type) i_iter->first)[*g_iter].second;
    683:             const int       size     = indSet.size();
    684:             int            *_indices = new int[size];
    685:
    686:             if (debug > 1) {std::cout << "  field " << *g_iter << std::endl;}
    687:             for(int i = 0; i < size; ++i) {
    688:               _indices[i] = indSet[i];
    689:               if (debug > 1) {std::cout << "    indices["<<i<<"] = " << _indices[i] << std::endl;}
    690:             }
    691:             disc->setIndices(_indices, i_iter->second);
    692:           }
    693:         }
    694:       };
    695:     public:
    696:       void setupField(const Obj<PETSC_MESH_TYPE> mesh, const Obj<PETSC_MESH_TYPE::real_section_type>& s, const int cellMarker = 2, const bool noUpdate = false, const bool setAll = false) {
    697:         typedef ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> Visitor;
    698:         const Obj<names_type>& discs  = this->getDiscretizations();
    699:         const int              debug  = s->debug();
    700:         names_type             bcLabels;
    701:
    702:         s->setChart(mesh->getSieve()->getChart());
    703:         int maxDof = this->setFiberDimensions(mesh, s, discs, bcLabels);
    704:         this->calculateIndices(mesh);
    705:         this->calculateIndicesExcluded(mesh, s, discs);
    706:         this->createReorder();
    707:         mesh->allocate(s);
    708:         s->defaultConstraintDof();
    709:         const Obj<PETSC_MESH_TYPE::label_type>& cellExclusion = mesh->getLabel("cellExclusion");
    710:
    711:         if (debug > 1) {std::cout << "Setting boundary values to " << std::endl;}
    712:         for(names_type::const_iterator n_iter = bcLabels.begin(); n_iter != bcLabels.end(); ++n_iter) {
    713:           Obj<PETSC_MESH_TYPE::label_sequence> boundaryCells = mesh->getLabelStratum(*n_iter, cellMarker);
    714:           if (setAll) boundaryCells = mesh->heightStratum(0);
    715:           //const Obj<PETSC_MESH_TYPE::real_section_type>&  coordinates   = mesh->getRealSection("coordinates");
    716:           const Obj<names_type>&         discs         = this->getDiscretizations();
    717:           const PETSC_MESH_TYPE::point_type               firstCell     = *boundaryCells->begin();
    718:           const int                      numFields     = discs->size();
    719:           PETSC_MESH_TYPE::real_section_type::value_type *values        = new PETSC_MESH_TYPE::real_section_type::value_type[mesh->sizeWithBC(s, firstCell)];
    720:           int                           *dofs          = new int[maxDof];
    721:           int                           *v             = new int[numFields];
    722:           Visitor pV((int) pow(mesh->getSieve()->getMaxConeSize(), mesh->depth())+1, true);
    723:
    724:           for(PETSC_MESH_TYPE::label_sequence::iterator c_iter = boundaryCells->begin(); c_iter != boundaryCells->end(); ++c_iter) {
    725:             ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mesh->getSieve(), *c_iter, pV);
    726:             const Visitor::point_type *oPoints = pV.getPoints();
    727:             const int                  oSize   = pV.getSize();
    728:
    729:             if (debug > 1) {std::cout << "  Boundary cell " << *c_iter << std::endl;}
    730:             //mesh->computeElementGeometry(coordinates, *c_iter, v0, J, NULL, detJ);
    731:             this->setCell(mesh, *c_iter);
    732:             for(int f = 0; f < numFields; ++f) v[f] = 0;
    733:             for(int cl = 0; cl < oSize; ++cl) {
    734:               //if (*c_iter == 0) std::cout << oPoints[cl] << std::endl;
    735:               const int cDim = s->getConstraintDimension(oPoints[cl]);
    736:               int       off  = 0;
    737:               int       f    = 0;
    738:               int       i    = -1;
    739:
    740:               if (debug > 1) {std::cout << "    point " << oPoints[cl] << std::endl;}
    741:               if (cDim || setAll) {
    742:                 if (debug > 1) {std::cout << "      constrained excMarker: " << mesh->getValue(cellExclusion, *c_iter) << std::endl;}
    743:                 for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
    744:                   const Obj<GeneralDiscretization>& disc    = this->getDiscretization(*f_iter);
    745:                   const Obj<names_type>           bcs     = disc->getBoundaryConditions();
    746:                   const int                       fDim    = s->getFiberDimension(oPoints[cl], f);//disc->getNumDof(this->depth(oPoints[cl]));
    747:                   const int                      *indices = disc->getIndices(mesh->getValue(cellExclusion, *c_iter));
    748:                   int                             b       = 0;
    749:
    750:                   for(names_type::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter, ++b) {
    751:                     const Obj<GeneralBoundaryCondition>& bc    = disc->getBoundaryCondition(*bc_iter);
    752:                     const int                            value = mesh->getValue(mesh->getLabel(bc->getLabelName()), oPoints[cl]);
    753:                     if (b > 0) v[f] -= fDim;
    754:                     if (value == bc->getMarker()) {
    755:                       if (debug > 1) {std::cout << "      field " << *f_iter << " marker " << value << std::endl;}
    756:                       for(int d = 0; d < fDim; ++d, ++v[f]) {
    757:                         dofs[++i] = off+d;
    758:                         //if (!noUpdate) values[indices[v[f]]] = (*bc->getDualIntegrator())(v0, J, v[f], bc->getFunction());
    759:                         if (!noUpdate) {
    760:                           values[indices[v[f]]] = bc->integrateDual(v[f]);
    761:                         }
    762:                         if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
    763:                       }
    764:                       // Allow only one condition per point
    765:                       ++b;
    766:                       break;
    767:                     } else {
    768:                       if (debug > 1) {std::cout << "      field " << *f_iter << std::endl;}
    769:                       for(int d = 0; d < fDim; ++d, ++v[f]) {
    770:                         if (!setAll) {
    771:                           values[indices[v[f]]] = 0.0;
    772:                         } else {
    773:                           //TODO: do strides of the rank of the unknown space; we need a natural way of handling vectors in a single discretization.
    774:                           values[indices[v[f]]] = bc->integrateDual(v[f]);
    775:                         }
    776:                         if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
    777:                       }
    778:                     }
    779:                   }
    780:                   if (b == 0) {
    781:                     if (debug > 1) {std::cout << "      field " << *f_iter << std::endl;}
    782:                     for(int d = 0; d < fDim; ++d, ++v[f]) {
    783:                       values[indices[v[f]]] = 0.0;
    784:                       if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
    785:                     }
    786:                   }
    787:                   off += fDim;
    788:                 }
    789:                 if (i != cDim-1) {throw ALE::Exception("Invalid constraint initialization");}
    790:                 s->setConstraintDof(oPoints[cl], dofs);
    791:               } else {
    792:                 if (debug > 1) {std::cout << "      unconstrained" << std::endl;}
    793:                 for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
    794:                   const Obj<GeneralDiscretization>& disc    = this->getDiscretization(*f_iter);
    795:                   const int                       fDim    = s->getFiberDimension(oPoints[cl], f);//disc->getNumDof(this->depth(oPoints[cl]));
    796:                   const int                      *indices = disc->getIndices(mesh->getValue(cellExclusion, *c_iter));
    797:
    798:                   if (debug > 1) {std::cout << "      field " << *f_iter << std::endl;}
    799:                   for(int d = 0; d < fDim; ++d, ++v[f]) {
    800:                     values[indices[v[f]]] = 0.0;
    801:                     if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
    802:                   }
    803:                 }
    804:               }
    805:             }
    806:             if (debug > 1) {
    807:               for(int f = 0; f < numFields; ++f) v[f] = 0;
    808:               for(int cl = 0; cl < oSize; ++cl) {
    809:                 int f = 0;
    810:                 for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
    811:                   const Obj<GeneralDiscretization>& disc    = this->getDiscretization(*f_iter);
    812:                   const int                       fDim    = s->getFiberDimension(oPoints[cl], f);
    813:                   const int                      *indices = disc->getIndices(mesh->getValue(cellExclusion, *c_iter));
    814:
    815:                   for(int d = 0; d < fDim; ++d, ++v[f]) {
    816:                     std::cout << "    "<<*f_iter<<"-value["<<indices[v[f]]<<"] " << values[indices[v[f]]] << std::endl;
    817:                   }
    818:                 }
    819:               }
    820:             }
    821:             if (!noUpdate) {
    822:               mesh->updateAll(s, *c_iter, values);
    823:             }
    824:             pV.clear();
    825:           }
    826:           delete [] dofs;
    827:           delete [] values;
    828:         }
    829:         if (debug > 1) {s->view("");}
    830:       };
    831:     };
    832:     /*

    836:     */

    841:     PetscErrorCode RHS_FEMProblem(::Mesh mesh, SectionReal X, SectionReal section,  void * ctx) {
    842:       GenericFormSubProblem * subproblem = (GenericFormSubProblem *)ctx;
    843:       Obj<PETSC_MESH_TYPE> m;

    847:       MeshGetMesh(mesh, m);
    848:       SectionRealZero(section);
    849:       Obj<PETSC_MESH_TYPE::real_section_type> s;
    850:       SectionRealGetSection(section, s);
    851:       // Loop over cells
    852:       //loop over integrals;

    854:       GenericFormSubProblem::names_type integral_names = subproblem->getIntegrals();
    855:       GenericFormSubProblem::names_type::const_iterator n_iter = integral_names.begin();
    856:       GenericFormSubProblem::names_type::const_iterator n_iter_end = integral_names.end();
    857:       //const Obj<PETSC_MESH_TYPE::order_type>& order = m->getFactory()->getGlobalOrder(m, "default", s);
    858:       while (n_iter != n_iter_end) {
    859:         Obj<GeneralIntegral> cur_integral = subproblem->getIntegral(*n_iter);
    860:         //get the integral's topological objects.
    861:         std::string cur_marker_name = cur_integral->getLabelName();
    862:         int cur_marker_num = cur_integral->getLabelMarker();
    863:         int cur_rank = cur_integral->getTensorRank();
    864:         int cur_dimension = cur_integral->getSpaceDimension();
    865:         Obj<PETSC_MESH_TYPE::label_sequence> integral_cells = m->getLabelStratum(cur_marker_name, cur_marker_num);
    866:         PETSC_MESH_TYPE::label_sequence::iterator ic_iter = integral_cells->begin();
    867:         PETSC_MESH_TYPE::label_sequence::iterator ic_iter_end = integral_cells->end();
    868:         if (cur_rank == 1) {
    869:           PetscScalar * values;
    870:           PetscMalloc(cur_dimension*sizeof(PetscScalar), &values);
    871:           PetscScalar * elemVec;
    872:           PetscMalloc(cur_dimension*sizeof(PetscScalar), &elemVec);
    873:           Obj<GenericFormSubProblem::names_type> discs = subproblem->getDiscretizations();
    874:           //loop over cells
    875:           while (ic_iter != ic_iter_end) {
    876:             subproblem->setCell(m, *ic_iter);
    877:             //loop over discretizations
    878:             GenericFormSubProblem::names_type::iterator d_iter = discs->begin();
    879:             GenericFormSubProblem::names_type::iterator d_iter_end = discs->end();
    880:             while (d_iter != d_iter_end) {
    881:               Obj<GeneralDiscretization> disc = subproblem->getDiscretization(*d_iter);
    882:               //evaluate the RHS at each index point
    883:               for (int i = 0; i < disc->size(); i++) {
    884:                 values[disc->getIndices()[i]] = -1.*disc->evaluateRHS(i);
    885:               }
    886:               d_iter++;
    887:             }
    888:             cur_integral->tabulateTensor(elemVec, values);
    889:             SectionRealUpdateAdd(section, *ic_iter, elemVec);
    890:             ic_iter++;
    891:           }
    892:           SectionRealGetSection(section, s);
    893:         } else if (cur_rank == 2) {
    894:           //cancel out BCs if necessary... this doesn't require any knowledge of the discretization form (if handled right).
    895:           PetscScalar * full_tensor;
    896:           PetscMalloc(cur_dimension*cur_dimension*sizeof(PetscScalar), &full_tensor);
    897:           PetscScalar * elemVec;
    898:           PetscMalloc(cur_dimension*sizeof(PetscScalar), &elemVec);
    899:           while (ic_iter != ic_iter_end) {
    900:             subproblem->setCell(m, *ic_iter);
    901:             cur_integral->tabulateTensor(full_tensor);
    902:             //create the linear contribution from the BCs
    903:             PetscScalar * xValues;
    904:             SectionRealRestrict(X, *ic_iter, &xValues); //get the coefficients -- BUG: setCell restricts as well; static
    905:             for(int f = 0; f < cur_dimension; f++) {
    906:               elemVec[f] = 0.;
    907:               for(int g = 0; g < cur_dimension; g++) {
    908:                 elemVec[f] += full_tensor[f*cur_dimension+g]*xValues[g];
    909:               }
    910:             }
    911:             SectionRealUpdateAdd(section, *ic_iter, elemVec);
    912:             ic_iter++;
    913:           }
    914:           //delete full_tensor;
    915:           //delete elemVec;
    916:         } else {
    917:           throw Exception("RHS_FEMProblem: Unsupported tensor rank");
    918:         }
    919:         n_iter++;
    920:         SectionRealGetSection(section, s);
    921:       }
    922:       // Exchange neighbors
    923:       SectionRealComplete(section);
    924:       // Subtract the constant
    925:       if (m->hasRealSection("constant")) {
    926:         const Obj<PETSC_MESH_TYPE::real_section_type>& constant = m->getRealSection("constant");
    927:         Obj<PETSC_MESH_TYPE::real_section_type>        s;
    928:
    929:         SectionRealGetSection(section, s);
    930:         s->axpy(-1.0, constant);
    931:       }
    932:       PetscBool  flag;
    933:       PetscOptionsHasName(NULL, "-vec_view", &flag);
    934:       if (flag) {
    935:         SectionRealGetSection(section, s);
    936:         s->view("RHS");
    937:       }
    938:       return(0);
    939:     }

    944:     PetscErrorCode Jac_FEMProblem(::Mesh mesh, SectionReal section, Mat A, void * ctx) {
    946:       GenericFormSubProblem * subproblem = (GenericFormSubProblem *)ctx;
    947:       Obj<PETSC_MESH_TYPE::real_section_type> s;
    948:       Obj<PETSC_MESH_TYPE> m;
    950:       MatZeroEntries(A);
    951:       MeshGetMesh(mesh, m);
    952:       SectionRealGetSection(section, s);
    953:       //loop over integrals; for now.
    954:       GenericFormSubProblem::names_type integral_names = subproblem->getIntegrals();
    955:       GenericFormSubProblem::names_type::iterator n_iter = integral_names.begin();
    956:       GenericFormSubProblem::names_type::iterator n_iter_end = integral_names.end();
    957:       const Obj<PETSC_MESH_TYPE::order_type>&        order         = m->getFactory()->getGlobalOrder(m, "default", s);
    958:       while (n_iter != n_iter_end) {
    959:         Obj<GeneralIntegral> cur_integral = subproblem->getIntegral(*n_iter);
    960:         //get the integral's topological objects.
    961:         std::string cur_marker_name = cur_integral->getLabelName();
    962:         int cur_marker_num = cur_integral->getLabelMarker();
    963:         int cur_rank = cur_integral->getTensorRank();
    964:         int cur_dimension = cur_integral->getSpaceDimension();
    965:         //GOOD LORD >:(
    966:         m->setMaxDof(subproblem->localSpaceDimension());

    968:         //divide here; ignore integral ranks that aren't 2 in the jacobian construction.
    969:         if (cur_rank == 2) {
    970:           Obj<PETSC_MESH_TYPE::label_sequence> integral_cells = m->getLabelStratum(cur_marker_name, cur_marker_num);
    971:           PETSC_MESH_TYPE::label_sequence::iterator ic_iter = integral_cells->begin();
    972:           PETSC_MESH_TYPE::label_sequence::iterator ic_iter_end = integral_cells->end();
    973:           PetscScalar * tensor;
    974:           PetscMalloc(cur_dimension*cur_dimension*sizeof(PetscScalar), &tensor);
    975:           while (ic_iter != ic_iter_end) {
    976:             subproblem->setCell(m, *ic_iter);
    977:             cur_integral->tabulateTensor(tensor);
    978:             updateOperator(A, m, s, order, *ic_iter, tensor, ADD_VALUES);
    979:             ic_iter++;
    980:           }
    981:         }
    982:         n_iter++;
    983:       }
    984:       MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    985:       MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
    986:       //MatView(A, PETSC_VIEWER_STDOUT_SELF);
    987:       return(0);
    988:     }

    993:     PetscErrorCode SubProblemView(SectionReal section, std::string name, PetscViewer viewer, int firstField = 0, int lastField = 0) {
    994:       //"vectorize" takes the first n discretizations and writes them out as a vector

    998:       Obj<PETSC_MESH_TYPE> m;
    999:       Obj<PETSC_MESH_TYPE::real_section_type> field;
    1000:       SectionRealGetBundle(section, m);
    1001:       SectionRealGetSection(section, field);
    1002:       const ALE::Obj<PETSC_MESH_TYPE::numbering_type>& numbering = m->getFactory()->getNumbering(m, 0);
    1003:       PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", numbering->getGlobalSize());

    1005:       if (lastField - firstField > 0) {
    1006:         PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name.c_str());
    1007:
    1008:       } else {
    1009:         if (name == "") {
    1010:           PetscViewerASCIIPrintf(viewer, "SCALARS Unknown double %d\n", 1);
    1011:         } else {
    1012:           PetscViewerASCIIPrintf(viewer, "SCALARS %s double %d\n", name.c_str(), 1);
    1013:         }
    1014:         PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
    1015:       }
    1016:       typedef PETSC_MESH_TYPE::real_section_type::value_type value_type;
    1017:       const Obj<PETSC_MESH_TYPE::real_section_type::chart_type>& chart   = field->getChart();
    1018:       const MPI_Datatype                  mpiType = ALE::New::ParallelFactory<value_type>::singleton(field->debug())->getMPIType();
    1019:       int enforceDim;
    1020:       int fiberDim = lastField - firstField + 1;
    1021:       if (lastField - firstField > 0) {
    1022:         enforceDim = 3;  //we need at least three vector components to be written out
    1023:       } else {
    1024:         enforceDim = 0;
    1025:       }
    1026:       if (field->commRank() == 0) {
    1027:         for(PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator p_iter = chart->begin(); p_iter != chart->end(); ++p_iter) {
    1028:           if (!numbering->hasPoint(*p_iter)) continue;
    1029:           const value_type *array = field->restrictPoint(*p_iter);
    1030:           const int&        dim   = field->getFiberDimension(*p_iter);
    1031:           ostringstream     line;
    1032:
    1033:           // Perhaps there should be a flag for excluding boundary values
    1034:           if (dim != 0) {
    1035:             for(int d = firstField; d <= lastField; d++) {
    1036:               if (d > 0) {
    1037:                 line << " ";
    1038:               }
    1039:               line << array[d];
    1040:             }
    1041:             for(int d = fiberDim; d < enforceDim; d++) {
    1042:               line << " 0.0";
    1043:             }
    1044:             line << std::endl;
    1045:             PetscViewerASCIIPrintf(viewer, "%s", line.str().c_str());
    1046:           }
    1047:         }
    1048:         for(int p = 1; p < field->commSize(); p++) {
    1049:           value_type *remoteValues;
    1050:           int         numLocalElementsAndFiberDim[2];
    1051:           int         size;
    1052:           MPI_Status  status;
    1053:
    1054:           MPI_Recv(numLocalElementsAndFiberDim, 2, MPI_INT, p, 1, field->comm(), &status);
    1055:           size = numLocalElementsAndFiberDim[0]*numLocalElementsAndFiberDim[1];
    1056:           PetscMalloc(size * sizeof(value_type), &remoteValues);
    1057:           MPI_Recv(remoteValues, size, mpiType, p, 1, field->comm(), &status);
    1058:
    1059:           for(int e = 0; e < numLocalElementsAndFiberDim[0]; e++) {
    1060:             for(int d = 0; d < fiberDim; d++) {
    1061:               if (mpiType == MPI_INT) {
    1062:                 PetscViewerASCIIPrintf(viewer, "%d", remoteValues[e*numLocalElementsAndFiberDim[1]+d]);
    1063:               } else {
    1064:                 PetscViewerASCIIPrintf(viewer, "%G", remoteValues[e*numLocalElementsAndFiberDim[1]+d]);
    1065:               }
    1066:             }
    1067:             for(int d = fiberDim; d < enforceDim; d++) {
    1068:               PetscViewerASCIIPrintf(viewer, " 0.0");
    1069:             }
    1070:             PetscViewerASCIIPrintf(viewer, "\n");
    1071:           }
    1072:           PetscFree(remoteValues);
    1073:         }
    1074:       } else {
    1075:         value_type *localValues;
    1076:         int         numLocalElements = numbering->getLocalSize();
    1077:         const int   size = numLocalElements*fiberDim;
    1078:         int         k = 0;
    1079:
    1080:         PetscMalloc(size * sizeof(value_type), &localValues);
    1081:         for(PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator p_iter = chart->begin(); p_iter != chart->end(); ++p_iter) {
    1082:           if (!numbering->hasPoint(*p_iter)) continue;
    1083:           if (numbering->isLocal(*p_iter)) {
    1084:             const value_type *array = field->restrictPoint(*p_iter);
    1085:             //const int&        dim   = field->getFiberDimension(*p_iter);
    1086:
    1087:             for(int i = firstField; i <= lastField; ++i) {
    1088:               localValues[k++] = array[i];
    1089:             }
    1090:           }
    1091:         }
    1092:         if (k != size) {
    1093:           SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of values to send for field, %d should be %d", k, size);
    1094:         }
    1095:         int numLocalElementsAndFiberDim[2] = {numLocalElements, fiberDim};
    1096:         MPI_Send(numLocalElementsAndFiberDim, 2, MPI_INT, 0, 1, field->comm());
    1097:         MPI_Send(localValues, size, mpiType, 0, 1, field->comm());
    1098:         PetscFree(localValues);
    1099:       }
    1100:       return(0);
    1101:     }
    1102:   }  //namespace Problem
    1103: } //namespace ALE

    1105: #endif