Actual source code: FEMProblem.hh
petsc-3.4.5 2014-06-29
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