1: #ifndef __UFCProblem
4: /*
5: UFC Problem class for sieve; allows for UFC-assembled problems to be "bundled"
7: */
10: #include <sieve/Mesh.hh>
11: #include <sieve/FEMProblem.hh>
12: #include <ufc.h>
14: namespace ALE {
15: namespace Problem {
17: class UFCHook { //initialize all the data from the UFC forms so we don't need to keep redoing it
18: public: 20: //the model for this is of course b(u, v) = (l(v), f)... because it's UFC.
21: //if we need to extend this make these into sets.
23: ufc::form * _lform;
24: ufc::form * _bform;
26: //finite elements
28: ufc::finite_element * _b_finite_element; //the thing upon which the cell integral is defined; however we want to crack it open for discretizations
29: ufc::finite_element * _l_finite_element; //ditto
31: std::map<int, ufc::finite_element *> _b_sub_finite_elements;
32: std::map<int, ufc::finite_element *> _l_sub_finite_elements;
34: ufc::shape _cell_shape;
36: ufc::cell * _cell; //assume, for now, that the mesh has one cell type
38: //cell integrals
40: std::map<int, ufc::cell_integral *> _b_cell_integrals;
41: std::map<int, ufc::cell_integral *> _l_cell_integrals;
43: //interior facet integrals
45: std::map<int, ufc::interior_facet_integral *> _b_interior_facet_integrals;
46: std::map<int, ufc::interior_facet_integral *> _l_interior_facet_integrals;
48: //exterior facet integrals
50: std::map<int, ufc::exterior_facet_integral *> _b_exterior_facet_integrals;
51: std::map<int, ufc::exterior_facet_integral *> _l_exterior_facet_integrals;
53: public: 55: UFCHook(ufc::form * bform, ufc::form * lform) {
57: _bform = bform;
58: _lform = lform;
59: 60: //set up the cell, the finite_elements, the cell integrals
61: 62: //create the cell integrals
63: //_num_l_cell_integrals = lform->num_cell_integrals();
64: //_l_cell_integrals = new ufc::cell_integral *[_num_l_cell_integrals];
65: for (unsigned int i = 0; i < lform->num_cell_integrals(); i++) {
66: _l_cell_integrals[i] = _lform->create_cell_integral(i);
67: }
69: if (_l_cell_integrals.size() == 0) throw Exception("No linear integrals in the form");
71: //_num_b_cell_integrals = bform->num_cell_integrals();
72: //_b_cell_integrals = new ufc::cell_integral *[_num_l_cell_integrals];
73: for (unsigned int i = 0; i < bform->num_cell_integrals(); i++) {
74: _b_cell_integrals[i] = bform->create_cell_integral(i);
75: }
76: 77: if (_b_cell_integrals.size() == 0) {
78: throw Exception("No bilinear integrals in the form");
79: }
80: 81: //create the finite elements for the bform and lform; they should be the same I do believe
83: _b_finite_element = _bform->create_finite_element(0); //the other appears to be for quadrature
85: //_num_finite_elements = _bform->num_coefficients() + _bform->rank();//_bform->num_finite_elements();
86: //_finite_elements = new ufc::finite_element *[_num_finite_elements];
88: if (_b_finite_element->num_sub_elements() > 1) for (unsigned int i = 0; i < _b_finite_element->num_sub_elements(); i++) {
89: _b_sub_finite_elements[i] = _b_finite_element->create_sub_element(i);
90: }
92: /*
93: if (_b_sub_finite_elements.size() == 0) {
94: throw Exception("No finite elements in the form");
95: }
96: */
98: _l_finite_element = _lform->create_finite_element(0);
100: //_num_coefficient_elements = _lform->num_coefficients() + _lform->rank();
101: //_coefficient_elements = new ufc::finite_element *[_num_coefficient_elements];
102: if (_l_finite_element->num_sub_elements() > 1) for (unsigned int i = 0; i < _l_finite_element->num_sub_elements(); i++) {
103: _l_sub_finite_elements[i] = lform->create_finite_element(i);
104: }
106: /*
107: if (_l_sub_finite_elements.size() == 0) {
108: throw Exception("No coefficient elements in the form");
109: }
110: */
112: int dim, embedDim;
113: //derive these from the mesh and/or the finite element;
114: _cell = new ufc::cell;
117: _cell_shape = _l_finite_element->cell_shape();
120: if (_cell_shape == ufc::interval) {
121: dim = 1;
122: embedDim = 1;
123: _cell->entity_indices = new unsigned int *[dim+1];
124: _cell->entity_indices[0] = new unsigned int[2];
125: _cell->entity_indices[1] = new unsigned int[1];
126: } else if (_cell_shape == ufc::triangle) {
127: dim = 2;
128: embedDim = 2;
129: _cell->entity_indices = new unsigned int *[dim+1];
130: _cell->entity_indices[0] = new unsigned int[3];
131: _cell->entity_indices[1] = new unsigned int[3];
132: _cell->entity_indices[2] = new unsigned int[1];
133: } else if (_cell_shape == ufc::tetrahedron) {
134: dim = 3;
135: embedDim = 3;
136: _cell->entity_indices = new unsigned int *[dim+1];
137: _cell->entity_indices[0] = new unsigned int[4];
138: _cell->entity_indices[1] = new unsigned int[6];
139: _cell->entity_indices[2] = new unsigned int[4];
140: _cell->entity_indices[3] = new unsigned int[1];
141: } else {
142: throw Exception("Unknown cell shape");
143: }
144: //create the cell
145: _cell->topological_dimension = dim;
146: _cell->geometric_dimension = embedDim;
147: double * tmpCoords = new double[(dim+1)*embedDim];
148: _cell->coordinates = new double *[dim+1];
149: for (int i = 0; i < dim+1; i++) {
150: _cell->coordinates[i] = &tmpCoords[i*embedDim];
151: }
152: 153: }
155: ~UFCHook() {
156: //remove all this stuff
157: delete _lform;
158: delete _bform;
159: std::map<int, ufc::cell_integral *>::iterator i_iter = _l_cell_integrals.begin();
160: std::map<int, ufc::cell_integral *>::iterator i_iter_end = _l_cell_integrals.end();
161: while (i_iter != i_iter_end) {
162: delete i_iter->second;
163: i_iter++;
164: }
165: i_iter = _b_cell_integrals.begin();
166: i_iter_end = _b_cell_integrals.end();
167: while (i_iter != i_iter_end) {
168: delete i_iter->second;
169: i_iter++;
170: }
171: std::map<int, ufc::finite_element *>::iterator f_iter = _l_sub_finite_elements.begin();
172: std::map<int, ufc::finite_element *>::iterator f_iter_end = _l_sub_finite_elements.end();
173: while (f_iter != f_iter_end) {
174: delete f_iter->second;
175: f_iter++;
176: }
177: delete _cell->coordinates[0];
178: delete _cell->coordinates;
179: }
181: //accessors
184: void setCell(Obj<PETSC_MESH_TYPE> m, PETSC_MESH_TYPE::point_type c) {
185: //setup the cell object such that it contains the mesh cell given
187: const Obj<PETSC_MESH_TYPE::sieve_type> s = m->getSieve();
188: Obj<PETSC_MESH_TYPE::real_section_type> coordinates = m->getRealSection("coordinates");
189: int dim = m->getDimension();
190: if ((int)_cell->topological_dimension != m->getDimension() - m->height(c)) throw ALE::Exception("Wrong element dimension for this UFC form");
191: //int depth = m->depth();
192: //the entity indices should be FIXED because we don't care to use any of the DoFMap stuff
193: ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV((int) pow(m->getSieve()->getMaxConeSize(), m->depth())+1, true);
194: ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*m->getSieve(), c, pV);
196: const PETSC_MESH_TYPE::point_type * oPoints = pV.getPoints();
197: int num_oPoints = pV.getSize();
198: int vertex_index = 0;
199: for (int t = 0; t < num_oPoints; t++) {
200: if (m->depth(oPoints[t]) == 0) {
201: // _cell->entity_indices[0][vertex_index] = oPoints[t];
202: const double * tmpcoords = coordinates->restrictPoint(oPoints[t]);
203: for (int i = 0; i < dim; i++) {
204: _cell->coordinates[vertex_index][i] = tmpcoords[i];
205: }
206: vertex_index++;
207: }
208: }
209: }
210: };
212: #if 0
214: class UFCCell : GeneralCell {
215: private:216: ufc::cell * _cell; //the UFC cell.
217: public:218: UFCCell (ufc::shape shape, int embed_dimension) : GeneralCell() {
219: _embedded_dimension = embed_dimension;
220: if (shape == ufc::interval) {
221: //allocate coordinates
222: _num_vertices = 2;
223: _coordinates = new double [_num_vertices*_embedded_dimension];
224: //setup the closure reorder
225: //setup the ufc::cell
226: } else if (shape == ufc::triangle) {
227: //allocate coordinates
228: _num_vertices = 3;
229: _coordinates = new double [_num_vertices*_embedded_dimension];
230: //setup the closure reorder
231: //setup the ufc::cell
232: } else if (shape == ufc::tetrahedron) {
233: //allocate coordinates
234: _num_vertices = 4;
235: _coordinates = new double [_num_vertices*_embedded_dimension];
236: //setup the closure reorder
237: //setup the ufc::cell
238: } else throw Exception("UFCCell: Unsupported Shape");
239: }
240: };
242: #endif
244: class UFCBoundaryCondition : public GeneralBoundaryCondition {
245: private:246: //exact function
248: ufc::function * _function; //this will be evaluated to set the conditions.
249: ufc::function * _exactSol; //this will be evaluated in the case that there is an analytic solution to the associated problem.
250: ufc::finite_element * _finite_element; //we have to evaluate both the function and the exactsol at some point
251: ufc::cell * _cell; //ugh
252: int * _closure2data; //map of LOCAL DISCRETIZATION closure dof indices to data dof indices; used for evaluation
255: public:256: UFCBoundaryCondition(MPI_Comm comm, int debug = 0) : GeneralBoundaryCondition(comm, debug) {
257: _function = NULL;
258: _exactSol = NULL;
259: }
261: UFCBoundaryCondition(MPI_Comm comm, ufc::function * function, ufc::function * exactsol, ufc::finite_element * element, ufc::cell * cell, const std::string& name, const int marker, int debug = 0) : GeneralBoundaryCondition(comm, debug){
262: _function = function;
263: _exactSol = exactsol;
264: this->setFiniteElement(element);
265: _cell = cell;
266: this->setLabelName(name);
267: this->setMarker(marker);
268: }
270: void setFiniteElement(ufc::finite_element * finite_element) {
271: _finite_element = finite_element;
272: }
274: ufc::finite_element * getFiniteElement() {
275: return _finite_element;
276: }
278: void setCell(ufc::cell * cell) {
279: _cell = cell;
280: }
282: ufc::cell * getCell() {
283: return _cell;
284: }
286: void setFunction(ufc::function * function) {
287: _function = function;
288: }
290: ufc::function * getFunction() {
291: return _function;
292: }
294: //SHOULD be the same as the closure reorder for the overlying discretization. This is just for library interface of course.
296: virtual void setReorder(int * reorder) {
297: this->_closure2data = reorder;
298: }
300: virtual const int * getReorder() {
301: return this->_closure2data;
302: }
305: //yeargh! do we ever use evaluate explicitly?
306: virtual double integrateDual(unsigned int dof) {
307: int ufc_dof;
308: ufc_dof = this->getReorder()[dof];
309: //just evaluate the function using the coordinates; assume the order doesn't get that screwed up
310: return this->_finite_element->evaluate_dof(ufc_dof, *_function, *_cell);
311: }
312: };
314: class UFCDiscretization : public GeneralDiscretization {
315: //implementation of the discretization class that provides UFC hooks.
316: private:317: //Each form will have the finite element -- not true; you have to associate by which part of the SPACE the
318: //integral exists over. So, we must assign each integral to the appropriate space (as we might have TWO integrals
319: //over say, the interior and then a boundary integral NOT TRUE AGAIN! all the integrals have to happen over the
320: //overall form; they must be done a class up. This should just be used to segment the space
321: //
322: ufc::cell * _cell;
323: ufc::finite_element * _finite_element; //this will typically be a subelement.
324: ufc::function * _rhs_function; //TODO: generalize this further
326: public:327: UFCDiscretization(MPI_Comm comm, int debug = 1) : GeneralDiscretization(comm, debug) {
328: }
330: UFCDiscretization(MPI_Comm comm, ufc::finite_element * element, ufc::cell * cell, int debug = 1) : GeneralDiscretization(comm, debug){
331: _finite_element = element;
332: _cell = cell;
333: }
335: ~UFCDiscretization() {
336: }
338: virtual void createReorder() {
340: //create the forward and backward maps to and from UFC data layout to closure data layout.
341: //do this case-by-case with info from the discretization
343: int * reorder = new int[this->size()];
344: std::map<int, int> closure_offsets;
345: int offset = 0;
346: if (this->getFiniteElement()->cell_shape() == ufc::interval) {
347: //find the indices corresponding to every point in the topology and assume they're ordered right; if not then we redo this later.
348: closure_offsets[0] = 0;
349: offset += this->getNumDof(1);
350: closure_offsets[1] = offset;
351: offset += this->getNumDof(0);
352: closure_offsets[2] = offset;
353: offset += this->getNumDof(0);
354: //we know the reorder; do it!
355: int offset = 0;
356: for (int i = 0; i < this->getNumDof(0); i++, offset++) reorder[closure_offsets[1] + i] = offset;
357: for (int i = 0; i < this->getNumDof(0); i++, offset++) reorder[closure_offsets[2] + i] = offset;
358: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[0] + i] = offset;
359: } else if (this->getFiniteElement()->cell_shape() == ufc::triangle) {
360: closure_offsets[0] = 0;
361: offset += this->getNumDof(2);
362: closure_offsets[1] = offset;
363: offset += this->getNumDof(1);
364: closure_offsets[2] = offset;
365: offset += this->getNumDof(1);
366: closure_offsets[3] = offset;
367: offset += this->getNumDof(1);
368: closure_offsets[4] = offset;
369: offset += this->getNumDof(0);
370: closure_offsets[5] = offset;
371: offset += this->getNumDof(0);
372: closure_offsets[6] = offset;
373: offset += this->getNumDof(0);
374: //we know the reorder; do it!
375: int offset = 0;
376: //this order is checked -- and appears to be right
377: for (int i = 0; i < this->getNumDof(0); i++, offset++) reorder[closure_offsets[4] + i] = offset;
378: for (int i = 0; i < this->getNumDof(0); i++, offset++) reorder[closure_offsets[5] + i] = offset;
379: for (int i = 0; i < this->getNumDof(0); i++, offset++) reorder[closure_offsets[6] + i] = offset;
380: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[2] + i] = offset;
381: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[3] + i] = offset;
382: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[1] + i] = offset;
383: for (int i = 0; i < this->getNumDof(2); i++, offset++) reorder[closure_offsets[0] + i] = offset;
384: //TET!
385: } else if (this->getFiniteElement()->cell_shape() == ufc::tetrahedron) {
386: closure_offsets[0] = 0;
387: offset += this->getNumDof(3);
388: closure_offsets[1] = offset;
389: offset += this->getNumDof(2);
390: closure_offsets[2] = offset;
391: offset += this->getNumDof(2);
392: closure_offsets[3] = offset;
393: offset += this->getNumDof(2);
394: closure_offsets[4] = offset;
395: offset += this->getNumDof(2);
396: closure_offsets[5] = offset;
397: offset += this->getNumDof(1);
398: closure_offsets[6] = offset;
399: offset += this->getNumDof(1);
400: closure_offsets[7] = offset;
401: offset += this->getNumDof(1);
402: closure_offsets[8] = offset;
403: offset += this->getNumDof(1);
404: closure_offsets[9] = offset;
405: offset += this->getNumDof(1);
406: closure_offsets[10] = offset;
407: offset += this->getNumDof(1);
408: closure_offsets[11] = offset;
409: offset += this->getNumDof(0);
410: closure_offsets[12] = offset;
411: offset += this->getNumDof(0);
412: closure_offsets[13] = offset;
413: offset += this->getNumDof(0);
414: closure_offsets[14] = offset;
415: offset += this->getNumDof(0);
416: //we know the reorder; do it!
417: int offset = 0;
418: for (int i = 0; i < this->getNumDof(0); i++, offset++) reorder[closure_offsets[11] + i] = offset;
419: for (int i = 0; i < this->getNumDof(0); i++, offset++) reorder[closure_offsets[12] + i] = offset;
420: for (int i = 0; i < this->getNumDof(0); i++, offset++) reorder[closure_offsets[13] + i] = offset;
421: for (int i = 0; i < this->getNumDof(0); i++, offset++) reorder[closure_offsets[14] + i] = offset;
422: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[10] + i] = offset;
423: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[9] + i] = offset;
424: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[6] + i] = offset;
425: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[8] + i] = offset;
426: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[7] + i] = offset;
427: for (int i = 0; i < this->getNumDof(1); i++, offset++) reorder[closure_offsets[5] + i] = offset;
428: for (int i = 0; i < this->getNumDof(2); i++, offset++) reorder[closure_offsets[4] + i] = offset;
429: for (int i = 0; i < this->getNumDof(2); i++, offset++) reorder[closure_offsets[3] + i] = offset;
430: for (int i = 0; i < this->getNumDof(2); i++, offset++) reorder[closure_offsets[2] + i] = offset;
431: for (int i = 0; i < this->getNumDof(2); i++, offset++) reorder[closure_offsets[1] + i] = offset;
432: for (int i = 0; i < this->getNumDof(3); i++, offset++) reorder[closure_offsets[0] + i] = offset;
433: } else {
434: throw Exception("UFCDiscretization->createReordering(): unsupported cell geometry");
435: }
436: this->setReorder(reorder);
437: //set this as the reorder of all boundary conditions involved here.
438: Obj<names_type> bcs = this->getBoundaryConditions();
439: names_type::iterator bc_iter = bcs->begin();
440: names_type::iterator bc_iter_end = bcs->end();
441: while (bc_iter != bc_iter_end) {
442: Obj<GeneralBoundaryCondition> bc = this->getBoundaryCondition(*bc_iter);
443: bc->setReorder(reorder);
444: bc_iter++;
445: }
446: }
448: virtual void setReorder(const int * _reordering) {
449: _closure2data = _reordering;
450: }
452: virtual const int * getReorder() {
453: return _closure2data;
454: }
456: void setRHSFunction(ufc::function * function) {
457: _rhs_function = function;
458: }
460: ufc::function * getRHSFunction() {
461: return _rhs_function;
462: }
464: virtual double evaluateRHS(int dof) {
466: //our notion of a degree of freedom in this case is a little different from the one used by UFC; they have the whole dimension of the vector space as
467: //one DOF with multiple coefficients; for the sake of being able to loop independent of the vector space here we split it up.
469: int ufc_dof;
470: if (_cell == NULL) throw Exception("UFCBoundaryCondition->integrateDual: cell is not initialized");
471: //TODO reordering
472: ufc_dof = this->getReorder()[dof];
473: //switch the order
475: //TODO: reorder based upon the cell attributes -- switch the dof_index to be the screwy indices as used by the
476: //cells... this requires knowledge of the dofs per dim here
478: //just evaluate the function using the coordinates; assume the order doesn't get that screwed up
480: return _finite_element->evaluate_dof(ufc_dof, *_rhs_function, *_cell);
481: }
483: int size() {
484: return _finite_element->space_dimension();
485: }
487: ufc::finite_element * getFiniteElement() {return _finite_element;};
488: void setFiniteElement(ufc::finite_element * element) {_finite_element = element;};
489: ufc::cell * getCell() {return _cell;};
490: void setCell(ufc::cell * cell) {_cell = cell;};
492: };
494: class UFCCellIntegral : public GeneralIntegral {
495: private:496: ufc::cell_integral * _integral;
497: ufc::cell * _cell;
498: const double ** _coefficients;
499: double * _tmp_tensor; //allows for us to reorder on the fly.
500: double * _tmp_coefficients;
501: public:502: UFCCellIntegral(MPI_Comm comm, int debug = 1) : GeneralIntegral(comm, debug) {
503: this->_tmp_tensor = NULL;
504: this->_tmp_coefficients = NULL;
505: this->_integral = NULL;
506: this->_cell = NULL;
507: this->_coefficients = NULL;
508: }
510: UFCCellIntegral(MPI_Comm comm, ufc::cell_integral * integral, ufc::cell * cell, name_type name, int label, int rank, int dimension, int ncoefficients = 0, int debug = 0) : GeneralIntegral(comm, name, label, ncoefficients, debug){
511: this->_integral = integral;
512: this->_cell = cell;
513: this->setTensorRank(rank);
514: this->setSpaceDimension(dimension);
515: this->setNumCoefficients(ncoefficients);
516: if (rank == 1) {
517: this->_tmp_tensor = new double[dimension];
518: } else if (rank == 2) {
519: this->_tmp_tensor = new double[dimension*dimension];
520: } else {
521: throw Exception("UFCCellIntegral::UFCCellIntegral(): unsupported tensor rank");
522: }
523: if (ncoefficients) {
524: this->_tmp_coefficients = new double[dimension*ncoefficients];
525: this->_coefficients = new const double*[ncoefficients];
526: } else this->_coefficients = NULL;
527: }
529: virtual ~UFCCellIntegral() {
530: delete _coefficients;
531: if (_tmp_tensor) delete _tmp_tensor;
532: if (_tmp_coefficients) delete _tmp_coefficients;
533: };
535: virtual ufc::cell * getCell() {
536: return _cell;
537: }
539: virtual void setCell(ufc::cell * cell) {
540: _cell = cell;
541: }
543: virtual void setIntegral(ufc::cell_integral * integral) {
544: _integral = integral;
545: }
547: virtual ufc::cell_integral * getIntegral () {
548: return _integral;
549: }
552: //there will eventually be some reordering occuring here, but for now just assume we've done the right thing.
553: virtual void tabulateTensor(double * tensor, const double * coefficients = NULL) {
554: //just run the frickin' tabulate_tensor routine from the form on what gets spit out by restrictClosure right now; rebuild or reorder later.
555: //the coefficients have come in in
556: if (this->getNumCoefficients()) {
557: for (int c = 0; c < this->getNumCoefficients(); c++) {
558: this->_coefficients[c] = &_tmp_coefficients[c*this->getNumCoefficients()];
559: for (int i = 0; i < this->getSpaceDimension(); i++) {
560: //reorder the coefficients into the _tmp_coefficients array
561: _tmp_coefficients[c*this->getNumCoefficients() + this->getReorder()[i]] = coefficients[i*this->getNumCoefficients() + c]; //may be wrong.
562: }
563: }
564: }
566: //the coefficients should already be reordered from the closure order to the normal order based upon the set
567: _integral->tabulate_tensor(_tmp_tensor, _coefficients, *_cell);
568: //reorder the tabulated tensor
569: if (this->getTensorRank() == 1) {
570: for (int f = 0; f < this->getSpaceDimension(); f++) {
571: tensor[f] = _tmp_tensor[this->getReorder()[f]];
572: }
573: }else if (this->getTensorRank() == 2) {
574: for (int f = 0; f < this->getSpaceDimension(); f++) {
575: for (int g = 0; g < this->getSpaceDimension(); g++) {
576: tensor[f*this->getSpaceDimension()+g] = _tmp_tensor[this->getReorder()[f]*this->getSpaceDimension()+this->getReorder()[g]]; //this should reorder it right.
577: }
578: }
579: }
580: }
581: };
583: //TODO: interior; exterior integral forms; these will have their own reordering notion.
585: class UFCFormSubProblem : public GenericFormSubProblem {
587: //Data and helper functions:
589: private:591: ufc::cell * _cell; //we'll have to pass this on to all the children of this class;
592: ufc::finite_element * _full_finite_element; //the big one with all unknowns (even mixed) described. all integrals are assumed to be defined on it
594: public:596: UFCFormSubProblem(MPI_Comm comm, int debug = 0) : GenericFormSubProblem(comm, debug) {
597: //yeah
598: _cell = NULL;
599: _full_finite_element = NULL;
600: }
602: UFCFormSubProblem(MPI_Comm comm, ufc::cell * cell, ufc::finite_element * finite_element, int debug = 0) : GenericFormSubProblem(comm, debug) {
603: _cell = cell; //yeargh!
604: _full_finite_element = finite_element;
605: }
607: UFCFormSubProblem(MPI_Comm comm, ufc::form * b_form, ufc::form * l_form, int debug = 0) : GenericFormSubProblem(comm, debug){
608: //set up the discretizations, integrals, and other things; discretization-level things like boundary conditions must be set separately of course.
609: //set the main finite element, the sub-elements, and
610: }
612: ufc::cell * getCell() {
613: return _cell;
614: }
616: ufc::finite_element * getFiniteElement() {
617: return _full_finite_element;
618: }
620: virtual int localSpaceDimension() {
621: return _full_finite_element->space_dimension();
622: }
624: virtual void createReorder () {
625: /*
626: GOAL: keep all interface-level things (in the other file) in *closure* order.
627: 629: The UFC ordering is assumed to go:
630: 1. discretization 1 entry 0 reordered per-cell like UFC
631: discretization 1 entry 1 ...
632: 2. discretization 2 reordered per-cell like UFC
633: ... etc
634: 635: So, we first pull apart each discretization and reorder (based upon the indices of that discretization and the rank of the form, then based upon the discretization
637: Then pass the overall reorder to each integral to make all of this go smoothly.
639: */
640: 641: int * reorder = new int[this->localSpaceDimension()];
642: //PetscPrintf(PETSC_COMM_WORLD, "starting the reorder calculation\n");
643: //loop over discretizations, creating the local reorders
644: Obj<names_type> discs = this->getDiscretizations();
645: int offset = 0;
646: for (names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); d_iter++) {
647: const Obj<GeneralDiscretization> disc = this->getDiscretization(*d_iter);
648: disc->createReorder();
649: const int * disc_indices = disc->getIndices();
650: const int * disc_reorder = disc->getReorder();
651: int disc_size = disc->size();
652: //PetscPrintf(PETSC_COMM_WORLD, "disc %s size: %d\n", (*d_iter).c_str(), disc_size);
653: for (int i = 0; i < disc_size; i++) {
654: reorder[disc_indices[i]] = disc_reorder[i] + offset;
655: //PetscPrintf(PETSC_COMM_WORLD, "map closure dof index %d to %d\n", disc_indices[i], reorder[disc_indices[i]]);
656: }
657: offset += disc_size;
658: //use the discretizations closure indices to reorder the whole-field indices
659: }
660: //set the reorder and set all the cell integrals to have this reorder as well. Unfortunately this will take some hacking for the interior facet integrals.
661: this->setReorder(reorder);
662: Obj<names_type> integrals = this->getIntegrals();
663: names_type::iterator i_iter = integrals->begin();
664: names_type::iterator i_iter_end = integrals->end();
665: while (i_iter != i_iter_end) {
666: Obj<GeneralIntegral> integral = this->getIntegral(*i_iter);
667: integral->setReorder(reorder);
668: i_iter++;
669: }
670: }
672: virtual void setCell(Obj<PETSC_MESH_TYPE> mesh, PETSC_MESH_TYPE::point_type c) const {
673: if (_cell == NULL) throw Exception("UFCFormSubProblem: Uninitialized cell");
674: Obj<PETSC_MESH_TYPE::real_section_type> coords = mesh->getRealSection("coordinates");
675: const double * cellcoords = mesh->restrictClosure(coords, c);
676: int index = 0;
677: for (unsigned int i = 0; i < _cell->topological_dimension+1; i++) {
678: for (unsigned int j = 0; j < _cell->geometric_dimension; j++) {
679: _cell->coordinates[i][j] = cellcoords[index];
680: index++;
681: }
682: }
683: }
685: ~UFCFormSubProblem(){};
687: };
688: }
689: }
691: #endif