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

      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